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ABSTRACT 

This paper investigates the hydrodynamic performances of an SPH code incorporating an artificial heat conductivity term in which 
the adopted signal velocity is applicable when gravity is present. To this end, we analyze results from simulations produced using a 
suite of standard hydrodynamical test problems. In accordance with previous findings it is shown that the performances of SPH to 
describe the development of Kelvin-Helmholtz instabilities depend strongly on the consistency of the initial condition set-up and on 
the leading error in the momentum equation due to incomplete kernel sampling. On the contrary, the presence of artificial conductivity 
does not significantly affect the results. 

An error and stability analysis shows that the quartic B-spline kernel (M 5 ) possesses very good stability properties and we propose 
its use with a large neighbor number, between ~ 50 (2D) to ~ 100 (3D), to improve convergence in simulation results without being 
affected by the so-called clumping instability. Moreover, the results of the Sod shock tube test demonstrate that in order to obtain 
simulation profiles in accord with the analytic solution, for simulations employing kernels with a non-zero first derivative at the origin 
it is necessary the use of a much larger number of neighbors than in the case of the M 5 runs. 

SPH simulations of the blob test show that in order to achieve blob disruption it is necessary the presence of an artificial conductivity 
term. However, it is found that in the regime of strong supersonic flows an appropriate limiting condition, which depends on the 
Prandtl number, must be imposed on the artificial conductivity SPH coefficients in order to avoid an unphysical amount of heat 
diffusion. 

Results from hydrodynamic simulations that include self-gravity show profiles of hydrodynamic variables that are in much better 
agreement with those produced using mesh-based codes. In particular, the final levels of core entropies in cosmological simulations 
of galaxy clusters are consistent with those found using AMR codes. This demonstrate that the proposed diffusion scheme is capable 
of mimicking the process of entropy mixing which is produced during structure formation because of diffusion due to turbulence. 
Finally, results of the Rayleigh-Taylor instability test demonstrate that in the regime of very subsonic flows the code has still several 
difficulties in the treatment of hydrodynamic instabilities. These problems being intrinsically due to the way in which in standard 
SPH gradients are calculated and not to the implementation of the artificial conductivity term. To overcome these difficulties several 
numerical schemes have been proposed which, if coupled with the SPH implementation presented in this paper, could solve the issues 
which recently have been addressed in investigating SPH performances to model subsonic turbulence. 
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1. INTRODUCTION 

Smoothed particle hydrodynamics (SPH) is a Lagrangian, mesh- 
free, particle method which is used to model fluid hydro- 
dynamics in numerical simulations. The techni que wa s orig- 
inally developed in an ast rophysical context dLucv I Il977t 
iGingold & Monaghanlll977 ), but since then i t has been widely 
applied in many other areas ( Monaghan I2005I) of computational 
fluid dynamics. 

The method has several properties which make its 
use particularly a dvantageous in astrophysical problem s 
(iHernquist & Katz I Il989t iRosswoel 120091 ISpringell l2010h . 
Because of its Lagrangian nature the development of large mat- 
ter concentrations in collapse problems is followed naturally, 
moreover the method is free of advection errors and is Galilean 
invariant, Finally, the method naturally incorporates self-gravity 
and possess very good conservation properties, its fluid equa- 
tions being derived from variational principles. 

The other computational fluid dynamical method com- 
monly employed in numerical astrophysics is the standard 
Eulerian grid based approach in which the fluid is evolved 



on a discretized mesh dStone & Norman I 1992t Rvu et al 



1993 



iNorman & Bryan I Il999bt iFrvxell et al. I l200ol: iTevssier 



2002). The spatial resolution of an Eulerian scheme based on 
a fixed Cartesian grid is often insufficient however to ade- 
quately resolve the large dynamic range frequently encoun- 
tered in many astrophysical problems, such as galaxy forma- 
tion. This has motivated the development of adaptative mesh 
refinement (AMR) methods, in which the spatial resolution of 
the grid is locally refined according to some selection criterion 
(iBerger & Colella 1 119891: Ikravtsov. Klvpin & Khokhlovlll997t 
INorman 112005 ). Additionally, the order of the numerical scheme 
has been i mproved by adopting the para bolic piecewise method 
(PPM) of IColella& Woodward I d 19841). such as in the AMR 
Eulerian codes ENZ O dNorman & Brvanlll999bl) and FLASH 
jFrvxeUetal.ll2000h . 

Application of these different types of hydrodynami- 
cal codes to the same test problem with identical ini- 
tial conditions should in principle lead to similar results. 
However, there has been growing evidence over recent 
years that in a variety of hydrodynamical test cases there 
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are significant difference s between the results produced the 
two t ypes of methods ([O'Shea et al. I l2005ak lAgertz et al. ' 



ypej 

2007; Wadsley, Veeravalli & Couch man 
2008t IMitchell et al. I l2009t |Read. Havfield & Agertz 



2008; 



Valcke et al. 11201 filJunk etal 



For instance, lAgertz et al 



120101). 



Tasker etak I 
2010L 



(l2007t) showed that in the stan- 



dard SPH formulations the growth of Kelvin-Helmholtz (KH) 
instabilities in shear flows is artificially suppressed when steep 
density gradients are present at the contact discontinuities. 
Moreover, the level of central entropy produced in binary 
merger non-radiative simulations of galaxy clusters is signifi- 
cantly higher by a factor ~ 2 in Eulerian simulatio ns than in 
those made using SPH codes ( IMitchell et al. 1 120091). The ori- 
gin of these discrepancies has been recognized ( Agertz et al. I 
120071: IRead. Havfield & Agertz l2010l) as partly due to the intrin- 
sic difficulty for SPH to properly model density gradients across 
fluid interfaces, which in turn implies the presence of a surface 
tension effect which inhibits the growth of KH instabilities. A 
second problem is due to the Lagrangian nature of SPH codes, 
which prevents mixing of flui d elements at the particle level and 
leads to entropy gen eration dWadslev. Veeravalli & Couchman I 
I2008t IMitchell et al. Il2009l) . In particular. IMitchell et al. I d2009t) 
argue that the main explanation for the different levels of central 
entropy found in cluster simulations is due to the different de- 
gree of entropy mixing present in the two codes. Unlike SPH, in 
Eulerian codes fluid mixing occurs by definition at the cell level 
and a certain degree of overmixing is certa inly present in th ose 
simulations made using mesh-based codes dSpringelll20i0bl) . 

Given the advantages of SPH codes highlighted pre- 
viously, it appears worth pursuing any improvement in 
the SPH method capable of overcoming its present limi- 
tations. Along this line of investigatio n many efforts have 
been 
2008 



made by a number of authors ( Abell 



2008 



2011: Price 



Wadsley, Veeravalli & Couchman 
20 lOt iRead. Havfield & Agertz I l20ldnHefi & Springe! 
ChaJjiutsuka & Navakshin 11201 OtlMurante et al. 11201 ll). 



Valcke et af 
2010t 



Price 1 (120081 ) introduces an artificial heat conduction term in 



the SPH equations with the purpose of smoothing thermal en- 
ergy at fluid interfaces. This artificial conductivity (AC) term 
in turn gives a smooth entropy transition at contact discon- 
tinuities with the effect of enforcing pressure continuity and 
removing the artificial surface tension effect which inhibits 
the growth of KH instabilities a t fluid interfaces. Similarly, 
IWadslev. Veeravalli & Couchman I d2008l) suggest that in SPH 
the lack of mixing at the particle level can be alleviated by 
adding a heat diffusion term to the equations so as to mimic the 
effects of subgrid turbulence, thereby improving the amount of 
mi xing. 

IRead. Havfield & Agertz] d2010l) present an SPH imple- 
mentation in which a m odified density estimate is adopted 
( Ritc hie & Thomas 1199 ll) . together with the use of a peaked ker- 
nel and a much larger number of neighbors. The authors showed 
that the new scheme is capable of following the develop ment of 
fluid i nstabilities in a better way than in standard SPH. lAbelll 
d201 lb presents an alternative derivation of the SPH force equa- 
tion which avoids the problem encountered by standard SPH in 
handling fluid instabilities, although the approach is not inher- 
ently energy or momentum conserving and is prone to large in- 
tegration errors when the simulati on resolutio n is low . 

The method proposed by I Inutsukal d2002l) reformu- 
lates the SPH equations by introducing a kernel convolu- 
tion so as to consistently calculate density and hydrodynamic 
forces. The latter are determined using a Riemann solver 
(Godunov SPH). The method has recently been revisited by 



ICha. Inutsuka & Navakshin I d2010l) and iMurante et al. I (1201 lh . 
who showed that the code correctly follows the development of 
fluid instabilities in a variety of hydrodynamic tests. 

A deeper modificat ion than those present ed here so far 
has been introduced by iHefi & Springel I d2010l) . who replaced 
the traditional SPH kernel approach with an new density esti- 
mate based on Voronoi tesselation. The authors showed that the 
method is free of surface tension effects and therefore the growth 
rate of fluid instabilities is not adversely affected as in standard 
SPH. 

Finall y, a radically new numerical scheme has been intro- 
duced by iSpringell d2010bl) . with the purpose of retaining the 
advantages of both SPH and mesh-based codes. In the new code 
the hydrodynamic equations are solved on a moving unstruc- 
tured mesh using a Godunov method with an exact Riemann 
solver. The mesh is defined by the Voronoi tesselation of a set 
of discrete points and is allowed to move freely with the fluid. 
The method is therefore adaptative in nature and thus Galilean 
invariant but, at the same time, the accuracy with which shocks 
and contact discontinu ities are described i s that of an Eulerian 
code. In a recent paper lBauer & Springell d2012l) argue that the 
standard formulation of SPH fails to accurately resolve t he de- 
velopm ent of turbulence in the subsonic regime (but see iPrice I 
d2012bl) for a different viewpoint). The authors draw their con- 
clusions by analyzing results from simulations of driven sub- 
sonic turbulence made using the new moving-mesh code, named 
AREPO, and a standard SPH code. 

S imilar conclusions were reached in a set of c ompanion pa- 
pers dSiiacki et al. 1201 lllVogelsberger et al. 1201 lh . in which the 
new code was used in galaxy formation studies to demonstrate 
its superiority over standard SPH. However, the code is charac- 
terized by considerable complexity which makes the use the SPH 
scheme still appealing and, more generally, it is desirable that 
simulation results produced with a specific code should be repro- 
duced with a completely independent numerical scheme when 
complex non-linear phenomena are involved. It appears worth- 
while, therefore, to investigate, along the line of previous au- 
thors, the possibility of constructing a numerical scheme based 
on the traditional SPH formulation which is capable of correctly 
describing the development of fluid instabilities and at the same 
time incorporates the effects of self-gravity when present. 

This is the aim of the present study, in which the SPH scheme 
is modified by incorp orating into th e equations an AC diffusion 
term as described by IPrice I (120081) . However, in IPrice I (120081) 
the strength of the AC is governed by a signal velocity which is 
based on pressure discontinuities. For simulations where grav- 
ity is present this approach is not applicable because hydro- 
static equilibrium requires pressure gradients. An appropriate 
signal velocity for conductivity when gravity is present is then 
used to construct an AC-SPH code with the purpose of treat- 
ing the growth of fluid instabilities self-consistently. The viabil- 
ity of the approach is tested using a suite of test problems in 
which results obtained using the new code are contrasted with 
the corresponding ones produced in previous work using differ- 
ent schejries ; _TjTe_£ode_2^ in form to that presented 
by [Wadslev, V eeravalli & Couchman I d2008l) . but here the en- 
ergy diffusion equation is implemented in a different manner. 

The paper is organized as follows. Sect. [2] presents the hy- 
drodynamical method and introduces the AC approach. In Sect. 
[3]we investigate the effectiveness of the method by presenting re- 
sults from a suite of purely hydrodynamic al test problems. These 
are: the two-dimensional Kelvin-Helmholtz instability, the Sod 
shock tube, the point explosion or Sedov blast wave test and the 
blob test. In Sect. [4] we then discuss results from hydrodynamic 
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tests that include self-gravity. Specifically, we consider the cold 
gas sphere or Evrard collapse test, the Rayleigh-Taylor instabil- 
ity and the cosmological integration of galaxy clusters. Finally, 
the main results are summarized in Sect. [5] 



2. The Hydrodynamic method 

Here we present the basic f eatures of the me thod , for recent re - 
views see lRosswogl d2009l) . ISpringel1 d2010h and lPricel d2012h . 



2.1. Basic SPH equations 

In the SPH method the fluid is described by a set of N par- 
ticles with mass m„ velocity v,, density p„ and a thermody- 
namic variable such as the specific thermal energy m, or the 
entropic function A,. The particle pressure is then defined as 
Pi = (y - l)p,M; = Aip y r where y = 5/3 for a monoatomic 
gas. The density estimate p(r) at the particle position r, is given 
by 



Pi = Y i mjW{\r i j\,hi), 



(1) 



where r (J = r, - rj, WQrij], hi) is the interp olating kerne l that 
has compact support and is zero for |ry| > (hi dPrice II20T2I) . The 
kernel is normalized to the condition J Wd D r = 1 . The sum in 
Eq. ([TJ is over a finite number of particles and the smoothing 
length hi is a variable that is implicitly defined by the equation 



f D {£,hi) D pi = N sp hmi , 



(2) 



where D is the number of spatial dimensions, = n, 4n/3 
for D = 2, 3 and N sp h is the number of neighboring particles 
within a radius £7z, . This equation can be rewritten by defining 
hi in units of the mean interparticle separation 



hi = nirndpif 10 , 



(3) 



so that N 2 ° h = tt(£t]) 2 and N™ h = 4n((?]) 3 /3. The smoothing 
length hi is determined by solving the non-linear equation ©; 
note that N sp h does not necessarily need to be an integer but can 
take arbitrary v alues if 77 is used as the fundamental parameter 
determining hi dPricell2012l) . A kernel commonly employed is 
the M4, or cubic spline, which is zero for £ > 2. In three dimen- 
sions typical choices for N sp h lie in the range N sp h ~ 33 - 64, 
which for the M4 kernel corresponds to 77 1.-1.25. The equa- 
tion of motion for the SP H particles can be derived fr om the 
Lagrangian of the system (ISpringel II20T0I; iPrice 1120121) and is 
given by 



dVi y-> 



1 ViWij(hi) + —1-ViWijQij) 



(4) 



where the coefficients Q, are defined as 



dh ^ dWikQii) 



dpi 



dh; 



(5) 



These terms are present in the momentum equation be- 
cause the smoothing length itself is implicitly a function of 
the particle coordinates through Eq. ©. 



2.2. Artificial viscosity 

In SPH, the momentum equation must be generalized by adding 
an appropriate viscous force, which is aimed at correctly treating 
the effects of shocks. An artificial viscosity (AV) term is then 
introduced with the purpose of dissipating local velocities and 
preventing particle interpenetration at the shock locations. The 
new term is given by 



(61 



where the term W/j = \{W{rjj, hi) + W(rij, hj)) is the sym- 
metrized ker nel and n, , is the AV tenso r. In the SPH entropy 
formulation (ISpringel & Hernquist l2002t) . it is the entropy func- 
tion per particle A, which is integrated and its time derivative is 
calculated as follows 



dAi 
dt 



ly- 1 

1 Pi 



^mjUijnj.ViWij^^) 



(7) 



where v,y = v, - v,-. For the AV tensor we adopt here the 
form proposed by dMonaghan 1 1 1 997h in analogy with Riemann 
solvers: 



Pij 



where 



vf / = a + Cj - 3p u 



(8) 



(9) 



r-,j < but zero 



is the signal velocity and p,j = Vjj ■ rij/\r/j\ if v,j 
otherwise, so that Yl/j is non-zero only for approaching particles 
Here scalar quantities with the subscripts i and j denote arith- 
metic averages, c, is the sound speed of particle ;, the parameter 
cti regulates the amount of AV, and is a controlling factor that 
reduces the stren gth of AV in th e presence of shear flows. The 
latter is given by dBalsara II 1995b 



fi = 



IV ■ v\i 



IV- v|i + ]Vxv|i ' 



(10) 



where (V • v), and ( V x v),- are the st andard SPH estimates 
for divergence and curl dMonaghan II20051) . For pure shear flows 
I V x v\i » |V ■ v\j so that the AV is strongly damped. 

Using the signal velocity the Courant condition on the 
timestep of particle ;' reads 



A& 0.3- 



c,|v*. v 



(11) 



In the standard SPH formulation the viscosity parameter a,- 
which controls the strength of the AV is giv en by a,; = const = 
ffo, with ffo = 1 being a common choice dMonaghan 1 120051) . 
This scheme has the disadvantage that it generates viscous dis- 
sipation in regions of the flow that are not undergoing shocks. In 
order to reduce spurious viscosity effects iMorris & Monaghan I 
dl997h proposed making the viscous coefficient a,- time depen- 
dent so that it can evolve in time under certain local conditions. 
The authors proposed an equation of the form 



dai_ _ ag - a,,,,-,, 
dt Tj 



+ Si, 



(12) 
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where S, -is a source term and t, regulates the decay of a, to 
the floor value a mm away from shocks. For the source term the 
following expression is adopted 



Si = S fi(-(V ■ v)i, 0)(a max - ttj) , 



(13) 



which is constructed in such a way that it increases in the 
presence of shocks. The prefactor Sq is unity for y = 5/3 and 
the damping factor is inserted to accou nt for the presence of 
vortici ty. The original form proposed by iMorris & Monaghan I 
dl997l) has been refined by the factor (a max — a,) dRosswog et al. I 
2000), which has the advantage of being more sensitive to 
shocks than the original formulation and of preventing a, 
from becoming higher th an a prescribed value a max . Recently, 
Cullen & Dehnenl d201 ll) pr esented an improved version of the 
Morris & Monaghan I d 19971) AV scheme, employing as shock 
indicator a switch based on the time derivative of the velocity 
divergence. 

The decay parameter t, is of the form 



Ti = 



hi 

Ci Id 



(14) 



where Id is a dimensionless parameter which controls 
the decay time sca le. In a number of test simulations, 
dRosswog et al. 1120001) found that appropriate values for the pa- 
rameters a max ,amin, and Id are 1.5,0.05, and 0.2, respectively. 
In principle, the effects of numerical viscosity in regions away 
from shocks can be reduced by setting /,/ to higher values than 
Id = 0.2. In practice, the minimum time necessary to propa- 
gate through the resolution length /?,- sets the upper limit Id = 1. 
However, the time evolution of the viscosity parameter can be af- 
fected if very short damping timescales are imposed. Neglecting 
variations in the coefficients, the solution to Eq. dT2l at times 
t > ti„ can be written as 



= q; + (ai(tm) - qd exp 



where 

Ti 



1+SiTi 

and qi is a modified source term 

_ &min "t" S iTj(Y max 



(15) 



(16) 



(17) 



From Eq. ( fT3T > it can be seen that 07(f) ^ a mnx in the strong 
shock regime S ,r; » 1 but this condition is not satisfied if 
SjTj i 1. Therefore for mild shocks this implies that, if very 
short decay time scales are imposed, the peak value of 07(f) at 
the shock front might be below the AV strength necessary to 
properly treat shocks. 

To overcome this difficulty a modification to dT3l has been 
adopted (Valdarnini 2011, herefater VI 1) which, when lj = 1, 
compensates the smaller values of S ,t, with respect to the small 
Id regimes. This is equivalent to considering a higher value for 
a fflM , so that in Eq. dT3l a max is substituted by a max — > %a max . 

The correction factor £ has been calibrated using the shock 
tube problem as reference and requiring a peak value of ~ 
0.6 - 0.7 for the viscosity parameter at the shock front, as in 
the = 0.2 case. The results indicate that a normalization of 
the form £ = (Zjj/0.2) ' 8 for / ( / > 0.2 satisfies these constraints. 
This normalization has been validated in a number of other test 
problems showing that it is able to produce a peak value of the 



Table 1. KH parameters for the simulations. From top to bottom: 
simulation label, Mach number M and x-velocity V] of the high- 
density layer, KH time scale tkh 



label run 


1 


2 


3 


4 


5 


M 


0.2 


0.4 


0.6 


0.8 


1.0 


h'vl 


0.26 


0.52 


0.77 


1.0 


1.3 




1.23 


0.56 


0.37 


0.28 


0.22 



viscosity parameter at the shock location which is independent 
of the chosen value of the decay parameter lj. 

The time-dependent AV formulation of SPH has been shown 
to be effective in reducing the damping of turbulence due to 
the e ffects of numerical viscosity in simulations of galaxy clus- 
ters dDolag et aL| 2005[ VI 1 ) Moreover, it has been used in 
a recent paper dPrice ]]2012bh to argue that the conclusion of 
iBauer & Springell d2012l) about the difficulty for SPH codes in 
properly modeling the development of subsonic turbulence is 
based just on using a SPH code in its standard AV formulation. In 
the following, unless otherwise stated, simulations will be per- 
formed using a time-dependent AV; this will be fully specified 
by the set of parameters {a min , a max , Id). 



2.3. Artificial conductivity 

As already outlined in the Introduction, different formulations 
have been proposed for overcoming the problems encountered 
by standard SPH in the treatment of fl uid discontin uities. Here 
we follow the approach suggested by iPrice I d2008l) . who pro- 
posed adding a dissipative term to the thermal energy equa- 
tion for smoothing the energy across contact discontinuities. The 
presence of these dissipative terms introduces a smoothing of 
entropy at fluid interfaces with the effect of removing pressure 
discontinuities. 

The m otivations ly i ng at the basis of this approach were dis- 
cussed by iMonaghan I d 1 9971) . who showed how the momentum 
and energy equations in SPH must contain a dissipative term re- 
lated to jumps in the variables across characteristics, in analogy 
with the corresponding Riemann solutions. 

The artificial conductivity (AC) term for the dissipation of 
energy takes the form 



I AC y P'J 



(18) 



where is the signal velocity, which does not need to be 
the same as that used in the momentum equation, en = rijlrij 
and af is the AC parameter which is of order unity Q Eq. ( fl~8b 
represents the SPH analogue of a diffusion equation of the form 



Df c V 2 Ui 



AC 

where 



V 2 u i = 2Y j m j 



Hi - Uj c, ■ VWjj 

Pj r U 



(19) 



(20) 



1 Note that th ere is a misprint in the sign of the corresponding Eq. 28 
of IPrice I f200^l 
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Table 2. Main simulation parameters of the KH tests. From top 
to bottom: simulation label, kernel function, neighbor number, 
signal velocity used in Eq. dT8l ). setup of initial HCP lattice: 
U=unsmoothed, S=smoothed. For all of the runs the AV param- 
eters are {a min ,a max ,l d } = {0.1, 1.5, 1). 



Simulations 


SPH 


RHO 


LIQ 


LP 


CRT 


M5 


Kernel 


CS 


CS 


LIQ 


LIQ 


CRT 


M5 


N s 


32 


32 


32 


32 


32 


50 




grav 


grav 


grav 


pres 


grav 


grav 


P 


U 


S 


S 


S 


S 


S 



is the SPH expression for the L aplacian dBrookshaw II 1985b 
and, in analogy with the analysis of lLodato & Price 1 (1201 Oh for 
defining an equivalent physical viscosity coefficient for the SPH 
numerical viscosity, we can define Df c as a numerical heat dif- 
fusion coefficient given by 



(21) 



ter af which must be constructed so that diffusion of thermal 
energy is kept under control and is damped away from disconti- 
nuities. This can be achieved by introducing a switch similar to 
that devised for AV. The dissipation parameter is then evolved 
according to 



daf 
dt 



(22) 



the meaning of the terms being similar to that of Eq. (fT2b . In 
particular, for the decay timescale rr = /i,/Cc, we set here C = 
0.2 and set to zero the floor value af ■ . For the source term S c 

mm i 

the following expression is used 



S C i=fchi 



\Jui + s 



("ma, - «f ) . 



(23) 



where fc is a dimensionless parameter of order unity and 
the parameter s avoids divergences as w,- — > 0. This equa- 
tion differs from the corr e sponding sou rce term proposed by 
IPrice &Monaghanl d2005l) . lPricel d2008l) by the factor af nax - af 
which has been inserted here in analogy with Eq. dT3l for in- 
troducing a stronger response of the switch in the presence of 
discontinuities. The choice of values for the parameters fc and 
a£„ r depends on the probl em under consideration, for example 
IPrice & Monaghanl (120051) proposed fc =0.1. Here a number of 
numerical experiments has shown that the best results in terms 
of response of the AC parameter af to the presence of discon- 
tinuities are obtained by setting fc — 1 and af naK = 1.5, which 
we assume henceforth as reference values. Moreover it is found 
that significant benefits in terms of sharpness of the AC param- 
eter profile across the discontinuity are obtained by inserting 
the af nax - af term. In principle, the choice of the derivative 
term used to dete ct discontinuities is arbitra ry, but in practice a 
second derivative dPrice & Monaghan 12005 ) term ensures better 
sensitivity to sharp disco ntinuities in thermal energy than a first 
derivative dPricell2005ab . 

An example of a signal velocity specifically designed to re- 
move pressure gradients at contact discontinuities is that origi- 



nally introduced by IPrice I d2008l) for pure hydrodynamical sim- 
ulations 



(24) 



and further refined by IValcke et al~l (120101) . The ability of the 
new AC formulation of SPH to follow the development of KH 
instabilities using this expressi on for the signal velocity has been 
verified in a numbe r of tests dPrice 1 120081: IValcke et aT~ll2010t 
iMerlin et al. I l2010h . However, the disadvantage of the signal 
speed d24b is that it cannot be applied when self-gravity is con- 
sidered because in that case a pressure gradient is present at hy- 
drostatic equilibrium. Imagine, for example, a self-gravitating 
gas sphere in hydrostatic equilibrium in which there is a negative 
radial temperature gradient, with the gas temperature decreasing 
when moving outward from the center of the sphere. An appli- 
cation of the SPH equations, with the AC term of Eq. d22l i now 
using the signal velocity (l24l . will lead to a heat flux from the 
inner to the outer regions and, in the long term, to the develop- 
ment of an unphysical isothermal temperature profile. A signal 
velocity which avoids this difficulty is given by 



An important issue concerns the choice of the AC parame- AC 



Vij(grav) = |(v ; - vj) 



(25) 



which corresponds to the f ormula tion proposed by 
IWadslev. Veeravalli & Couchman I (120081) in which a dissi- 
pative term is added to the evolution of the thermal equation 
with the purpose of modeling heat diffusion due to turbulence. 
The new term is constructed i n analogy with the subgrid-scale 
model of ISmagorinskv I d 19631) and is given by 



5t- c ? 



nij\Vi -Vj\(hj + hj) 
PU 



(ui-Uj)ejj ■ V,^ 



(26) 



where C is a coefficient of order unity, whose pre- 
cise value depends on the problem under considera- 
tion. For the rising hot bubb le prob lem considered by 
IWadslev. Veeravalli & Couchman I ((2008), the best agree- 
ment with the corresponding PPM results taken as reference 
is obtained setting C = 0.1, higher values being too much 
diffusive. Throughout this paper the SPH simulations of the 
hydrodynamic tests are performed by adopting the expression 
(fTSl for the ther mal energy dissipative term. With r e spect to the 
formulation of IWadslev. Veeravalli & Couchman I d2008l) this 
approach presents the advantage of using a diffusion parameter 
which is not constant but evolves in time according to a source 
term, so that the amount of diffusion away from discontinuities 
is minimized. For the AC signal velocity, we then use expression 
(1251) . whose performances in the AC formulation dT8l has not 
been fully tested before in SPH simulations of hydrodynamic 
test problems. 

Note that, using the signal velocity d25l ). the Von Neumann 
stability criterion becomes unimportant with respect the Courant 
condition dTTT ). The Von Neumann constraint requires At < 
Q.5Ax z /D, which in SPH reads 



2D AC 



(27) 



Since a c ^ 0(1), this condition implies Atf c > Atf in both the 
supersonic and subsonic regimes. 
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Fig. 1. Density maps for some of the 2D KH instability tests described in Sect. 13. II From left to right, each column shows the panels 
for runs having initial conditions with the same Mach number at the corresponding time t = tkh, as listed in Table Q] From top to 
bottom, different rows refer to the following simulations : SPH, RHO, LIQ and LP. The latter use the linear quartic kernel but with 
the si gnal velocity d24|). w hereas the first three use the expression J25b (see Table [2]). The plots can be compared directly with Fig. 
10 of lValckeetaLlOOlOh . 



3. Hydrodynamic tests 

In the following, simulation results obtained by applying the 
new AC-SPH code to a number of hydrodynamic test prob- 
lems are discussed with the objective of validating the code 
and assessing its performance. To this end, the simulation re- 
sults of the tests will be compared with the corresponding ones 
obtained by previous authors using different codes and/or nu- 
merical methods. The problems considered are usually pre- 
sented in the literature in dimensional or complexity order, 
but here we follow a different approach. Since t he KH insta- 
bility is the hydrodynamic test which initially dAgertz et al. I 
120071) originated the discussion about the difficulties of stan- 
dard SPH in properly handling the development of instabil- 
ities, and it is also the most widely considered in this con- 
text (lAbelll 1201 It iPricel l2008t ICha. Inutsuka & Navakshinl 



20Tot IHeB & Springelll2010t iRead. Havfield & Agertz 1 120 lOt 
Valcke et al. 2010t IMurante et al. 1 201 lh . we here discuss first 



in detail the two-dimensional version of the test. The setup of 
the other runs will then be considered in the light of the results 
obtained from the 2D KH test. 



3.1. 2D Kelvin-Helmholtz instability 



The KH instability arises in the presence of shear flow between 
two fluid layers when a small velocity perturbation is imposed 
in the direction perpendicular to the interface between the two 
fluids. The development of the instability is characterized by an 
initial phase, where the fluids interpenetrate each other, and then 
the forming of vortices, which become progressively more pro- 
nounced and turn into KH rolls at the onset of non-linearity. 
In the case of incompressible fluids for a sinusoidal perturba- 
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tion of wavelength A, this phase is reached with a time-scale 
riChandrasekharll 19611) 



TKH 



J-(p\ + P2) 
(PIPl) 



1/2, 



(28) 



where p\ and po are the fluid densities and v = Vi - V2 is 
the relative shear velocity. A proper modeling of instabilities in 
numerical hydrodynamics is essential since the KH instability 
plays a crucial role in the development of turbulence which oc- 
curs in many hydrodynamic al phenomena. In particular, the KH 
instability is relevant in many astrophysical processes, such as 
gas stripping from cold gas clouds occurring in galactic haloes 
and the production of entropy in the intracluster gas of galaxy 
clusters due to injection of turbulence. 

The growth of the KH instability in hydrod ynamic 
simul atio ns has b een a ddr essed by various authors dAbell 
201 U iPrice I 120081: IWadslev. Veeravalli & Couchman 



Hefi &Springell 1201 Ot iCha. Inutsuka & Navakshinl I201C ; 
Murante et al. 11201 ll) . These studies found that the development 



2008t iRead. Havfield & Agertzl l2010t IValcke et al. I I 201C : 



of the instability is artificially suppressed because of two 
distinct effects: the first problem is the so-called local mixing 
instability (LMI) which is due to entropy conservation and 
which inhibits mixing at the kernel scale thus introducing 
pressure discontinuities; the second problem arises because of 
errors in the momentum equation which cannot be reduced by 
increasing the number of neighbor particles without causing 
particle clumping. 

Given the wide variety of numerical parameters and initial 
conditions with which the KH instability has been addressed, we 
choose here to pe rform the tests using as reference the simula- 
tions presented bv lValcke et al. I (1201 C ). In particular, the authors 
implement in their SPH equations an AC term as that of Eq. ( TT81 ). 
but with the AC parameter set to unity, and with a signal velocity 
given by Eq. (124b . The authors performed a systematic analysis 
of the capabilities of SPH to capture the KH instability using 
different SPH formulations, kernels, numerical resolutions and 
KH time-scales. This choice allows to assess code capabilities 
by constructing a large suite of simulations whose numerical re- 
su lts can be contrasted against the corresponding ones discussed 
bv lValcke et all (l2010h . 

3.1.1. Initial conditions set-up 

The problem domain consists of a periodic box of unit length 
with cartesian coordinates x e {0, 1), y € {0, 1). Within the do- 
main there is a fluid with adiabatic index y — 5/3 which satisfies 
the following conditions: 



p, T, v x = 



Pi,T u vi \y-Q.5\ < 0.25 
p2,T 2 ,v 2 \y - 0.5| > 0.25. 



(29) 



As in IValcke et al. I (l20Toh . we choose here p\ = 10 and 
P2 — 1, with the index 1 referring to the high density layer. 
This choice is motivated by the finding that the difficulties of 
SPH in reproducing KH instabilities increase as the density con- 
trast between the two fluid layers gets higher. The two layers are 
in pressure equilibrium with Pi = P2 = 10 so that the sound 
velocities in the two layers are c\ = yjyP\]f>\ = 1-29 and 
C2 = yJyPi/pi = 4.08 , respectively. The layers slide against 
each other with opposite shearing velocities v\ = —vi and in 



order to seed the KH instability a small single-mode velocity 
perturbation is imposed along the y— direction 



wo sin(27rjc//l) . 



(30) 



where wq = 0.025 and A — 1/6. In order to restrict the per- 
turbation to spatial regions in the proximity of the interfaces, the 
perturbation ( l30l l is applied only if \y - cr\ < 0.025, where cr 
takes the values of 0.25 and 0.75, respectively. With this choice 
of parameters, the Mach number of the high-density layer is 
M =a \>i/ci ^ vi/1.29 and the KH time-scale is t kh 0.29/vi. 
As in IValcke et al.l d2010l) . we consider KH simulations with a 
range of five different Mach numbers; Table[T]reports the values 
of M together with the respective values of vi and tkh- 

In order to implement the the density set-up d29l ), a two- 
dimensional lattice of equal mass particles is placed inside the 
simulation box. We adopt here an isotropic hexagonal-close- 
packed (HCP) configuration for the particles coordinates instead 
of the more commonly employed Cartesian grid. The advantage 
of this configuration is that, for a given number of neighbors, it 
gives a better density estimate due to its symmetry properties. To 
construct the initial density configuration d29l . the lattice spac- 
ing of the particles is varied until the SPH density estimate ([D 
satisfies the required values within a certain tolerance criterion 
( i 1%) for the relative density error. The simulations were run 
using a total number of N = 5 1 2 2 particles and the initial particle 
specific energies were assigned after the density calculation so 
as to satisfy pressure equilib rium. This particle n umber is larger 
than that used in the runs of IValcke et al. I d2010t) (- 2QQK), but 
guarantees a dens ity setup with the sp ecified tolerance criterion. 

As noticed bv lValcke et al. I (1201 Oh and other authors, SPH is 
a numerical method which can only represent smoothed quanti- 
ties, and so applying it to hydrodynamic problems where strong 
density gradients are present can lead to inconsistencies. This is, 
in fact, the situation encountered by standard SPH in the treat- 
ment of KH instabilities, where the lack of entropy mixing in- 
duces an artificial pressure discontinuity at fluid interfaces with 
a jump in density. 

Motivated by these difficulties, we consider here a set of runs 
in which the density discontinuity at the interfaces is replaced by 
a smooth transition. To allow making a cons i stent c omparison 
with the corresponding runs of IValcke et al. I d2010l) . we adopt 
the same smoothing profile 



p(y) = D ± Aatan [B(y' + Q] , 

where the coefficients are given by 

A = Pi-Pi 
2atan(jS) ' 

6 



C 



D 



_6 
~2 ' 

Pi +P2 



(31) 



(32) 



(33) 



(34) 



(35) 



and y = y - cr + 6/2. In Eq. (f3TT > the sign in front of the coef- 
ficient A refers to cr = 0.25, 0.75 respectively. The parameters B 
and 5 determine the thickness of the density transition and take 
the values B = 10, 6 = 0.5/0.75. 

The following procedure is adopted in order to construct a 
lattice configuration which satisfies the density profile given by 



7 



R. Valdarnini: Hydrodynamic of an SPH code incorporating an artificial conductivity term 









■ 




1 1 1 1 


I I 1 1 






■ 


- 








U D.4 0,9 M 1 



Fig. 2. As in Fig.Q]but for the set of simulations LIQ, CRT and M 5 . 



Eq. (fJTJi. An HCP lattice of particle is first constructed in the 
domain < x < 1 and < y < 0.5 with a spacing ad own 
such that the SPH density is p = p 2 - Proceeding upward from 
y — 0, the lattice spacing is progressively reduced according to 
a = adownlpilpiy)] 112 until y = 0.25 and A^ OH ,„ particles are left. 
The same procedure is applied to a high-density lattice which in 
the same domain satisfies the condition p = p\\ starting from 
y = 0.5 and proceeding downward, the lattice spacing is in- 
creased so that a = a up [pi/p(y)] 1/2 until y = 0.25 and N up par- 
ticles remain of the original high-density lattice. The whole pro- 
cedure is numerically iterated until the numbers Ndown an d N up 
satisfy the conditions N doml +N up -N/2 and N up /N down = p\ /p 2 . 
The lattice is then replicated in the top half of the box (y — > 1 -y) 
so that the initial conditions are fully symmetric around y = 0.5. 
In the following, simulations with initial conditions for which 
the particle positions are arranged in a unsmoothed HCP lattice 
are denoted with the suffix SPH, whereas all of the others adopt 
the smoothed density profile constructed according to the proce- 
dure described here. 

Another issue which is found to have a significant impact 
on the ability of standard SPH to properly capture KH in- 
stabilities is the choice of the kernel function. According to 
iRead. Havfield & Agertz I d2010l) . the accuracy of the momen- 
tum equation for particle ; is governed by the leading error 



p>_ + p± 

Pj pi 



hiViWij 



(36) 



which vanishes in the continuum limit. However, this does 
not hold for a finite number of irregularly distributed particles 



or, more specifically, in the proximity of a contact discontinu- 
ity where a density step is present. A possible solution con- 
sists of increasing the number of neighbors so as to improve 
kernel sampling, although this approach presents some difficul- 
ties when the commonly employed M4 or cubic spline (CS) 
SPH kernel is used. This occurs because for a large number of 
neighbors the CS kernel is subject to the so-called clumping in- 
stability, in which pairs of particles with interparticle distance 
q = r/h < 2/3 remain close together because the CS kernel 
gradient tends to zero below this threshold distance. A s t ability 
analys i s dMorrislll996t lB0rve. Oman g & Trulsen 1 12004c iPrice I 
l2005at IRead. Havfield & Agertz 11201 Oh showslhat for the CS 
kernel, a Cartesian lattice of particles is unstable for 77 ^ 1.5 
or N sp h ~ 28, 110 when D - 2 and D = 3, respectively. The 
clumping degrades the spatial resolution because it reduces the 
effective number of neighbors with which integrals are sampled, 
thus still having large Eq errors even when the resolution is in- 
creased. To overcome this problem one can modify the kernel 
shape in order to have a non-zero kernel derivative at the origin. 
However, for a fixed number of neighbors, kernels of this kind 
have the drawback of giving a less accurate density estimate than 
that given by the CS kernel since near the origin the kernels have 
a steeper profile. 

As an alternative to the CS kerne l, we consider here th e linear 
quartic (LIQ) kernel, introduced by IValcke et al. I d2010h . which 
is a quartic polynomial for q s < q < g — 2 and is linear for 
< q < q s . The choice of the parameter q s determines the qual- 
ity of density estima tes. From a set of 2D Sod shock tube runs, 
IValcke et al. I d2010l) recommend q s = 0.6, which is the value 
adopted h ere. For the function al form and normalization of the 
kernel see lValcke et al. I d2010h . 
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Fig. 3. Rendered plots of the a c parameters are shown for the same runs as in Fig. [2] The distribution of particle values af has been 
interpolated at the map grid points according to the SPH prescription. 



Another kernel which has been introduced for the purpose 
of a voiding particle clumping is the core triangle (CRT) ker- 
nel ( iRead. Havfield & Agertzl 1201 Oh which has constant first 
derivative for < q < a and is similar to the CS kernel for 
a < q < £ = 2. This kernel is of the form 



w(q) = h 



(-3a + |a 2 ) q + 1 + |a 2 - |a 3 < q < a; 
1 - |<7 2 + 4<7 3 a < q < 1; 

1 < q < 2; 
q > 2, 



h 2 + W 
k(2-qf, 
0. 



(37) 



where l/cr = 2n(± + £ - ±a 5 ), An(\ + £ - f£) for D = 
2 and D = 3 , respectively. The value of a is fixed by the con- 
dition of continuity for the second derivative, giving a = 2/3. 
For the grid of initial conditions previously described, the sam- 
ple of KH simulations is then constructed by performing SPH 
runs with the same initial conditions but using different ker- 
nels. We consider the CS kernel, together with the LIQ and 
CRT kernels. For all of these runs the number of neighbors is 
N sp h = 32 (t?-1.5). 

However, another solution for reducing sampling errors con- 
sists of s till keeping a B-spline kernel function but increasing 
its order dPrice 1120121 sect.5.4). This approach presents the ad- 
vantage of retaining the bell shape of the kernel, a feature which 



provides good density estimates dFulk & Ouinn 1 19961) . After the 
M\ (CS) kernel, at the next order is the Ms or quartic kernel 



w(q) 



cr 
h 5 



o. 



^-5(1 
qy-5(l-q) 



q) + 



10 (I -qf 



0<q< 
\<q< 



1 



> 



<q< 

5 
2- 

has 



(38) 



which is truncated to zero for £ > 2.5 and 
96/1 199tt, 1 /20tt for D = 2 and D = 3, respectively. 

An additional set of SPH simulations was then performed 
using the M5 spline as the kernel. For consistency with the other 
runs, we kept the same ratio of smoothing length to particle spac- 
ing (77 ~ 1.5) so that for these runs the chosen number of neigh- 
bors is N sp h = 50. A non trivial issue concerns the role of pairing 
instability for this class of kernels. Because the gradient of the 
M5 kernel still goes to zero as q — > 0, one would expect the in- 
stability still to occur for rj ~ 1.5. Nonetheless, it will be seen 
that KH simulations with the M5 kernel do not exhibit particle 
clumping, in contrast with corresponding simulations performed 
with the CS kernel. This suggests that the stability properties of 
the Ms kernel are better than those of the CS one; this topic will 
be discussed in a subsequent dedicated section (13.1.3b . 

Table [2] summarizes the main simulation parameters used in 
the KH tests. For comparison purposes, a set of mirror runs was 
carried out for the LIQ simulations in which Eq. (l25l l, for the 
signal velocity, was replaced in Eq. ( fTSl by Eq. (l24l i which is 
based on pressure discontinuities. 
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Fig. 4. As in Figures Q] and [2]but here at the time t 
with high Mach number. 



2 for all of the panels. Note from Table Q]that this implies t » tkh for runs 



3.1.2. Results 

Fig. [1] shows the density maps at t — tkh for some of the KH 
simulations listed in Table [2] The panels can be compared di- 
rectly with the corresponding ones in Fig. 10 of IValcke et al. I 
(l2010h . A visual inspection reveals that the expected features of 
the KH runs are correctly reproduced and of a general consis- 
tency between the results produced by the two codes. In partic- 
ular, the LP simulations which use the pressure-based AC signal 
velocities d24l i appear to produce results almost identical to the 
corresponding LIQ ones. Thus showing how, for the KH tests 
examined here, the use of the two signal velocities can be con- 
sidered equivalent, within the spanned range of Mach numbers. 

Note, however, the absence of KH rolls for th e LP run with 
M — 0.4, whereas for the same simulation in IValcke et al. I 
d2010h the rolls have already been developed. Given the general 
agreement between the two codes it is hard to ascertain the ori- 
gin of the discrepancy for this specific run. However, the panels 
of Fig.rjjshow that the AC-SPH formulation, and this point will 
be discussed in detail later in sect. 14.21 is still unable to properly 



capture the development of KH instabilities for very subsonic 
shear flows. This suggests how small differences in the time in- 
tegration procedure of the two codes can become manifest in the 
long-term evolution of cold flows. 

In Fig. |2] the results for the same set of KH simulations of 
Table \T\ are shown, but with the use of kernels LIQ, CRT and 
M5. An important feature is now the appearance for the CRT 
and Ms runs of KH rolls in the M = 0.4 case. This improvement 
in the description of KH instabilities suggests that integration 
errors, that are present with the LIQ kernel, are now absent or 
reduced in the CRT and Ms runs. However, in the M = 0.2 case, 
the kernels are still unable to capture the development of the 
KH instability. For this test case, a high-resolution run (N = 
1024 2 ) using the Ms kernel (not shown here) still reveals the 
absence of any KH roll at t — tkh, thus showing that in SPH 
the problem of an accurate description of KH instabilities in the 
very subsonic regimes is not a resolution issue or due to th e AC 
implementation. See also iMcNallv. Lyra & Passv I (1201 ll) for a 
recent discussion on this topic. 
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To test the effectiveness of the switch dT3l in ensuring a fast 
response of the af parameters to the presence of thermal energy 
discontinuities, Fig.[3]renders the plots of the AC parameters that 
are shown for the test runs of Fig. [2] The maps are color coded 
according to the range of values of a.. The highest values lie in 
the range ~ 0.7 - 0.8 and are confined in a narrow strip around 
the interface layers, with the floor value ~ as the background 
value away from the discontinuities. Note that in general, and in 
particular for the M = 1 test case, the maximum values of the a c 
for the Ms runs are below those of the other runs. A result due 
to a faster growth rate of the instability, ensured by an improved 
kernel accuracy in interpolating the variables. 

The long term evolution of the KH tests is shown at t — 2 
in Fig. |U The overall features of the density plots for differ- 
e nt runs are similar to their counterparts displayed in Fig. 11 
of lValcke et al. I d2010l) . Note the tendency for the M 5 runs to the 
appearance of small scale features at the layer contacts. 

The aim of this section is to test the consistency of the 
present AC implementation by comparing results, extracted 
from a suite of AC-SPH simulatio ns of the 2D KH in stability 
problem, with those of similar runs dValcke et al. Il2010h . The re- 
sults of the KH tests indicate a code behavior whic h is in accord 
with w hat expected and with the simulations of IValcke et al. I 

(EoToh . 

However. lValcke et al. I d2010h argue that it is not the absence 
of the AC term the main reason of the SPH failure to develop 
KH instabilities, although the presence of AC is necessary for 
the long-term evolution of the instabilities, but this difficulty of 
SPH is rather due to two distinct reasons. The first issue is a 
general problem of consistency of the initial condition set-up, as 
SPH by definition can only deal with smoothed quantities and 
therefore its application to problems where sharp discontinuities 
are present le ads to inconsistencies. T his is part of the more gen- 
eral problem dRobertson et al. I [20101) of properly smoothing in 
numerical simulations of hydrodynamic test problems the dis- 
continuities initially present at the interfaces, so to achieve con- 
vergence in the solution. This aspect of the KH test problem can 
be cured by properly stretching the initial particle lattice so as to 
introduce a smooth density transition at the interfaces. The re- 
sults indicate a significant improvement in the capability of SPH 
to properly capture the correct growth rate of the KH instability. 

The other issue which causes a poor performance of SPH 
when handling KH instabilities is the leading error in the mo- 
mentum equation due to incomplete kernel sampling. This er- 
ror can be reduced by removing particle clumping, which de- 
pends on the kernel stability pr operties. In fact, the kernels LIQ 
and CRT have been introduced dRead. Havfield & Agertzll2010t 
IValckeet al. l2010h with the aim of removing the clumping insta- 
bility. Since the results of the simulations suggest performances 
for the Ms kernel that are quite similar to those achieved by these 
kernels, it is therefore interesting to better quantify the relative 
performances of these kernels. 

Following lJunk et aT71 d2010h . we then measure the growth 
rate of the KH instability for some of the runs and compare 
it with the linear theory growth rate expectation oc e t < rKH . This 
is achieved, using Fourier transforms, by measuring the time- 
evolution of the A - 1/6 growing mode of the v v velocity per- 
turbation compon ent. The details of the procedure are given in 
lJunk et all d2010h and will not be reported here. For the sake of 
clarity, the results of the LIQ runs are not displayed since they 
are quite similar to those of the CRT simulations. Moreover, we 
only display the growth rates for three different KH test cases, 
those with Mach number M = 0.2, 0.4 and M = 1, the results of 
the others being intermediate between these. 



There are a number of distinct features which are appar- 
ent from the panels of Fig. [5] The first is that simulations with 
smoothed IC ( RHO) perform systematically better than the un- 
smoothed ones (SPH). The second is that only for M = 1 the 
simulations with kernels CRT and Ms are able to correctly re- 
cover the expected growth rate. Finally, this capability degrades 
progressively as one considers lower Mach numbers. While this 
behavior is in accordance with the visual impression of the maps 
previously displayed, and confirms that the present SPH imple- 
mentation still has problems in the description of instabilities in 
subsonic flows, it is interesting to note that the performances of 
the Ms runs are similar to those obtained using the CRT kernel. 
This behavior is particularly interesting, since the simulations 
have been performed setting the ratio of the smoothing lengths 
to particle spacing to r\ 1 .5 so that the pairing instability, that 
is present in the runs which use the CS kernel, should be present 
in the Ms simulations as well. 

How the clumping instability affects sampling errors can 
be assessed by computing the particle errors E9, according to 
Eq. d36b . The distribution, throughout the simulation domain of 
the E9 errors versus y, is shown in the top panels of Fig. [6] for 
those M — 1 simulations performed using different kernels. The 
solid lines represent the mea n of the binned distri butions. Similar 
plots have been produced bv lValcke et al. I d2010t) and their Fig. 3 
can be used for comparative purposes. 

As expected, the largest E° errors are present in the SPH 
simulation, whereas better results are obtained by using the 
LIQ and CRT kernels, as indicated by the error distribution in 
their respective panels. This is a result of the absence of parti- 
cle clumping for these simulati ons, owing to the specific stabil- 
ity properties of th ese kernels dRead. Havfield & Agerfz |[20 1 Qt 
IValcke et al.ll2010l) . 

A striking feature is given vice versa by the E° error distribu- 
tion of the simulation performed using the Ms kernel, which, in 
fact, shows how the magnitude of the E° errors are even smaller 
than those of the two runs CRT and LIQ for this kernel. A result 
which is in accordance with what has been found previously by 
analyzing the growth rates, suggesting that, since all of the simu- 
lations were performed using the same rj, the stability properties 
of the Ms kernel are better than those of the CS one. 

To further investigate this point, the bottom panels of Fig. [6] 
show the 'nearest neighbor' map of the corresponding top pan- 
els. This is defined by interpolating the quantity qj at the grid 
points, according to the SPH prescription, where for the k - th 
neighbor q k . = \x/ - Xk\/h, and the neighbors have been sorted 
according to their distance from the particle i itself. 

The map of the Ms kernel indicates a distribution of the sec- 
ond nearest neighbor distribution qj that is quite similar to that 
displayed by the CRT and LIQ kernels. The only difference be- 
ing at the layer interfaces where the distribution of quantities qj 
for the Ms kernel is slightly shifted towards smaller values than 
for the other kernels. Note, vice versa, how for the CS kernel 
because of the pairing instability, the distribution of the neigh- 
boring distances is flipped with respect that of the other kernels. 

The results of Fig. |6l however, clearly show the absence of 
clumping instability for the Ms kernel. In the next section, we 
investigate, in more detail, the stability properties of this kernel. 



3.1.3. Stability issues 

The stability properties of SPH have be e n inv e stigated 
by a number of author s dSwegle & Hicks I Il995t iMorris I 



1996h iMonaghanl |2000t lB0rve. Omang & Trulsen I |200 



iRead. Havfield & Agertz I I20 lot iDehnen&Alvl l2012h . The 
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Fig. 5. Time evolution of the velocity field amplitude in the y direction, as measured by the kA = I2n mode of the Fourier transform 
of v y , for some of the KH instability tests described in Sect. I3.ll Each panel refers to KH simulations performed with the same 
Mach number, the initial conditions set-up being given in TableQ] Within each panel, the different curves are for AC-SPH runs with 
different simulation parameters, as specified in Table|2] The black solid line indicates the expected linear theory growth rate. 




Fig. 6. Top panels: distribution at time t = tkh of the errors | plotted versus y, as defined by Eq. d36l l. for the KH runs with 
Mach number M — 1 and different simulation parameters (See Table [2]). The red histograms show the mean values of the binned 
distributions. Bottom panels: each panel shows the nearest neighbor map of the run in the corresponding top panel. This is defined 
as the distribution, interpolated at the map grid points of the normalized distances q 2 = 6? I hi, where 5* is the distance \x, - Xu\ of 
the k-th neighbor of the particle i and the neighbors are sorted so that <S? < . 



instabilities are studied analytically by analyzing the dispersion mass m, density po, pressure Po and sound speed c 2 = yPo/po is 
relation for sound waves of small amplitude, propagating in a perturbed by a small wave 
uniform medium. We now derive the dispersion relation for the 
SPH equation of motion in a manner similar to the analysis of 
previous authors. A uniform lattice of equal mass particles of 

Xi = Xi° + a exp[k ■ Xi° - cjt] , (39) 
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Fig. 7. Contour plots of the dispersion relation d40b are shown as a function of wavenumber k and smoothing length h for a regular 
lattice of particles with spacing A, we consider a wavevector with orientation k = £(1,0, 0). Gray areas indicate the instability 
regions for which a> 2 < 0. From the left to right, the top panels show in the case of the Ms kernel the instability regions of the 
longitudinal and the two transverse waves. The bottom panels refer to the CS (M4) kernel. 



where a is the perturbation, k the wavevector and xp are the un- 
perturbed particle positions. By linearizing the equation of mo- 
tion for the perturbation one has the dispersion relation 



(jj 2 a u 



2tnPo 



Pi 



^ H^il - cos(k ■ Xij°)a v + 



Po 



where the summations are over the neighbors j of particle i, 
the vector q- t is defined as 



qt = ^ sin A: ■ *;/V, W, 7 , 



and H(W) is the Hessian of the kernel 

H = d 2 W(r) 
/JV dx u dx v 



(40) 



(41) 



For a given smoothing length h and wavevector k, the sta- 
ble solutions of Eq. (l40l > are defined by the condition ur > . 
Solutions for which or < represent exponentially growing or 
decaying perturbations. Moreover, it is useful to define a numer- 
ical sound speed C 2 m = or Ik 2 and a scaled numerical sound 
speed c 2 mm = C 2 um /c s . Clearly, the condition c 2 mm = 1 should be 
satisfied to correctly model the sound speed. 

We now examine the stability properties of the CS and Ms 
kernels. For simplicity, we consider a plane wave propagating 
along the jc-axis, k = £(1,0,0), and assume y = 5/3. The bot- 
tom panels of Fig|7]show, for the longitudinal and the two trans- 
verse waves of the perturbation, the instability regions of the CS 
kernel, these are denoted by the gray areas and represent the so- 
lutions to Eq. ( |40T > , in the domain (h,k), for which or < 0. 
Similarly, the top panels are for the Ms kernel. 



The longitudinal wave perturbation is responsible for the 
clumping instability, whereas the traverse waves give rise to th e 
so-called banding instability dRead. Havfield & AgertzlfeOlOl) . 
Unlike the clumping instability , the latter is relatively unimpor - 
tant in causing sampling errors dRead. Havfield & Agertz 11201 Oh 
and will not be further considered here. A comparison of the two 
stability plots for the longitudinal wave solution clearly shows 
that the stability properties of the Ms kernel are much better than 
those of the CS one. 

This i s part of a mor e general result which was already recog- 
nized by dMorrislll996l) : the stability properties of SPH are im- 
proved, as higher order spline kernels are used in place of the CS 
kernel. This occurs, basically, because the higher the order of the 
spline, the better it approximates a Gaussian kernel. In Eq. d40b 
one can see that the numerical sound speed c 2 um depends on the 
first and second derivative of kernel W. Ideally, one should have 
c num - 1 to accurately describe the sound wave propagation and 
this condition is fulfilled as smoother kernels are used, so as to 
keep the numerical dependency of c 2 um as weak as possible. 

Note however that here, differently from the stability prop- 
erties of the CRT kernel, the clumping instability is not entirely 
suppressed but rather it is present whenever r\ > 2.5. 

Finally, it must be stressed that the paring instability that oc- 
curs in the hydrodynamic tests described here, is unlikely to have 
a significant impact on many astr ophysical problems of inter est, 
where very cold flows are absent. iRasio & Lombardil d 19991) es- 
timate for instance, from SPH simulations of a stationary fluid, 
that lattice effects become important when velocity dispersions 
are below ~ 3% - 4% the sound speed. 

The results of this section and of the previous one, therefore, 
suggest that in order to avoid the pairing instability, the Ms ker- 
nel can be considered a viable alternative to the use of the other- 
wise steeper CRT and LIQ kernels, provided that the parameter 
77 is consistently rescaled. In the next section, we then discuss 
the related performances of these kernels in a test simulation. 
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3.2. Sod's shock tube 

A classic test used to investigate the hydrod ynamic capabilities 
of SPH codes is the Sod shock tube problem dHernauist & Katz ^ 



19891: IWadslev. Stadel & Ouinnl 120041: ISpringel I l2005t [Price 



2008t iTasker et al. 112008b iRosswog Il2009h . This test consists in 
a fluid, initially at rest, in which a membrane located at x — 
separates the fluid on the right, with high density and pressure, 
from the fluid on the left, with lower density and pressure. The 
membrane is removed at t — and a shock wave develops prop- 
agating toward the left, followed by a contact discontinuity and 
a rarefaction wave propagating to the right. 

A well known problem of standard SPH codes in reproduc- 
ing the analytic solution of the shock tube problem is the pres- 
ence of a pressure discontinuity that arises across the propagat- 
ing contact discontinuity. Simulations incorporatin g an artificial 
thermal conductivity term in the SPH equations dPrice I 120081 
lRosswog"ll2009t iPrice II20T2I) show shock profiles in which den- 
sity and thermal energy are resolved across the discontinuity, 
hence giving a continuous pressure profile. However, in these 
runs the AC formulation adopts the AC signal velocity (l24l 
where jumps in thermal energy are smoothed in accordance with 
the presence of pressure discontinuities. This is in contrast with 
the AC signal velocity d25T l employed here, for which in absence 
of shear flows, contact discontinuities are unaffected and there- 
fore cannot be used to remove the blip seen at the contact dis- 
continuity in the pressure profile of the shock tube SPH runs. 

However, the application of the AC-SPH code to this test is 
nonetheless of interest because it can still be used to validate 
code performances. In particular, we consider a 3D setup of the 
shock tube test and with these initial conditions we construct a 
set of AC-SPH simulations performed with different kernels and 
neighbor numbers. Shock tube profiles extracted from these sim- 
ulations are compared with the aim of assessing the goodness of 
different kernels in reproducing, for a given test problem, pro- 
files of known analytic solutions. 

The initial condition setup consists of an ideal fluid with 
y = 5/3, initially at rest at t = 0. An interface set at the 
origin separates the fluid on the right with density and pres- 
sure (pi,P\) = (4, 1) from the fluid on the left with {p2,Pi) = 
(1,0.1795). To construct these initial conditions, a cubic box of 
side-length unity was filled with 10 6 equal mass particles, so that 
800, 000 were placed in the right-half of the cube and 200, 000 in 
the left-half. The particles were extracted from two independent 
uniform glass-like distributions contained in a unit box consist- 
ing of 1.6 ■ 10 6 and 4 • 10 5 particles, respectively. This initial 
condition setup is the same previously implemented in Sect. 5.1 
of VI 1, to test the time-dependent AV scheme described in Sect. 

to 



For the same initial setup, to investigate the performances of 
different kernels in reproducing the analytic profiles of the shock 
tube problem, we perform runs with different kernels and neigh- 
bor numbers. The kernels considered are : CS, LIQ, CRT and 
M5. The number of neighbors varies in power of two between 
32 and 128. The SPH runs were realized by imposing periodic 
boundary conditions along the y and z axes and the results were 
examined at t — 0.12. We show results obtained using the stan- 
dard AV formulation, the results of the other runs being unim- 
portant from the viewpoint of kernel performances. 

Some of the simulation profiles extracted from the 3D SPH 
runs of the shock tube test are shown in Fig. [8] A striking feature 
is the wide scatter between the pressure profiles of the runs per- 
formed using different kernels or neighbor numbers. The same 



behavior is present for the thermal energy profiles, whilst very 
similar density profiles are exhibited by the same runs. 

There are several conclusions that can be drawn from the 
results of Fig. [8] The performances of the run are quite 

similar to those of CS-32, although for the former simulation a 
closer inspection reveals a slightly better treatment of the ther- 
mal energy spike at the contact discontinuity, the spike being due 
to the initial condition set-up. 

The worst results are obtained by the LIQ simulations when 
using N sp h = 32 or N sp i, = 64 neighbors. The profiles of the 
LIQ- 128 run are not shown here to avoid overcrowding in the 
plots, and are quite similar to those of simulation CRT- 128. This 
clearly demonstrates the need for this class of kernels to use a 
large ( say i 128 ) number of neighbors, so as to compensate for 
the density underestimate due to the steeper profiles introduced 
to avoid the pairing instability. 

To better quantify this density bias, Table |3]reports, for dif- 
ferent kernels and neighbor number, the mean SPH density es- 
timated from a glass-like configuration of N = 10 6 particles, in 
a unit periodic box of total mass one. The results illustrate how 
the density estimate of the M5 kernel with 64 neighbors is com- 
parable to the CS one obtained using 32 neighbors, the value of 
r\ being the same (77 1). Vice versa, for the LIQ and CRT ker- 
nels, only when N sp i, = 128 does the mean density approach the 
M 5 -64 estimate. 

The density values given in Table [3] also help to explain 
the rather poor performances of the LIQ kernel when using 32 
neighbors. From Fig. [8] it can be seen that for the correspond- 
ing run the relative pressure error is sp 10% in the unper- 
turbed right zone of the cube. This error is already present in 
the pressure profile at t = and it is originated from the den- 
sity error due to the adopted kernel and neighbor number, to- 
gether with the use of an entropy-conserving code to perform 
the simulations. Initially, particle entropies are assigned accord- 
ing to the initial conditions so that, for particles satisfying x, > 0, 
Aj = A\ - P\lp\ — 0.1. During the integration particle pressures 
are calculated according to P, = A;p?, and for x > a relative 
pressure error sp - ys p - 10% is present when e p 6%. 

Taken at face value, the results of Table [3] demonstrate that 
in 3D simulations a conservative lower limit for the kernels 
with a modified shape, should be to a ssume at least N xn h ~ 200 
neighbors. In fact, in a recent paper iRead & Havfield I (1201 2|) 
presented a new formulation of SPH where they adopt ed, as 
referen ce, the so-called high-order core triangle (HOCT, IPrice I 
(12005 ah ). The profile of this kernel is a generalization of that of 
the CRT kernel and the authors assume N sp i, = 442 as the refer- 
ence value for their scheme. 

However, the results of the previous sections suggest that the 
stability properties of the M5 kernel can be profitably used, with 
an appropriate choice of the rj parameter, to avoid the pairing 
instability when dealing with tests of hydrodynamic instabilities 
in which cold flows are present. 

Finally, the density estimates of Table [3] suggest that great 
care should be used when deciding on the goodness of a partic- 
ular kernel on the basis of its relative performances in terms of 
density estimates. The results of the 3D Sod shock tube clearly 
indicate how there could be a large difference between the sim- 
ulation and the expected solution profile of some hydrodynamic 
variable, such as pressure or thermal energy, while having, at 
the same time, a much smaller difference in the corresponding 
density profile. 
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Fig. 8. Results at t — 0.12 of the 3D shock tube test. The profiles of: density, pressure, thermal energy and velocity are plotted 
clockwise from top left. The solid black line represents the analytical solution, while lines with different styles and colors are the 
profiles of the AC-SPH runs with different kernels and neighbor numbers, as illustrated in the pressure panel. 



3.3. Sedov blast wave 

The Sedov blast wave test is used to validate, in three dimen- 
sions, the code capability in the strong shock regime. The test 
consists of a certain amount of energy E which is injected at 
t — into a very small volume of an ambient medium of uniform 
density p. The spherically symmetric shock propagates outward 
from the initi al volume and at the time t the shock front is located 
at the radius (ISedov II 19591) 

R(t)^/3(Et 2 /p) 1/5 , (42) 

where /3 ~ 1.15 fory = 5/3. 

Previous inv estigations (Rosswog, & Price I 120071: 
iMerlin et al. I 1201 Ol) showed that, owing to the large dis- 
continuities initially present in the thermal energy, incorporating 
an artificial conduction term in the energy equation greatly 
improves the description of the shock front in the simula- 
tions. Without the presence of this term, the initially large 
discontinuity in the thermal energy will soon give rise to a 
disordered particle distrib ution thus degrading the shock profile 
dRosswog & Price Il2007h . 

The initial setting of the test is realized as follows. A HCP 
lattice of N - 2 ■ 64 3 equal mass particles is arranged in a cu- 
bic box of side length unity. The particle masses are chosen so 
as to give p — 1 and periodic boundary conditions are imposed. 
The nearest particle to the position {0.5, 0.5, 0.5} is chosen as the 



Table 3. Average densities and sample standard deviations esti- 
mated from the SPH densities of a configuration of one million 
particles. These are arranged in a glass-like configuration inside 
a cubic periodic box of side length unity and total mass one. 
The SPH particle densities have been calculated for a variety of 
kernels and neighbor numbers N sp h (see text). 







N sph 




Kernel 


32 


64 


128 


CS 


1.005 ± 0.035 


1.002 ± 0.019 


1.0018 ± 0.012 


LIQ 


1.06 ±0.028 


1.024 ±0.016 


1.01 ± 0.011 


CRT 


1.079 ±0.039 


1.035 ±0.020 


1.016 ±0.012 


M5 


1.022 ±0.052 


1.003 ±0.024 


1.0008 ±0.015 



center particle. In order to consistently represent a point-like ex- 
plosion with the given numerical resolution, the particles j com- 
prised within the kernel radius of the center particle i are 
given an initial thermal energy such that the total injected energy 
is E = 1 . This blast wave energy is distributed among the neigh- 
boring particles not uniformly but with a weight proportional to 
W U . 

Using this initial condition set-up, we ran three test sim- 
ulations, which differ in the choice of the adopted kernel and 
neighbor number. For the CS kernel, we use N sp h = 64 neigh- 
bors. We also run two other test cases, now using the CRT 
and M5 kernel and N sp h = 128 neighbors. All of the simula- 
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tions were performed using the implemented AC scheme and 
the time-dependent AV formulation with parameters given by 
{a min ,a max ,l d ) = {0.1, 1.5,0.2} . 

Fig. [9] shows the radial density profiles at t — 0.06 for the 
different test runs. The solid black line represents the expected 
analytic solution, with the shock front being located at r - 0.37. 
Radial simulation profiles were obtained by averaging, for each 
radial bin r*, SPH densities calculated from the particle distri- 
butions over a set of (6,(p) = (20,20) grid points uniformly 
spaced in angular coordinates: these are located at the surface 
of a sphere with radius r^. The radial spacing is not uniform but 
is chosen so as to guarantee an accurate sampling of density in 
the proximity of the shock front, with about ~ 40 radial bins 
between r ~ 0.34 and r ~ 0.42. 

All of the simulations are in fair agreement with the ana- 
lytical solution and the differences between the simulation pro- 
files are negligible. At the shock front, the profiles exhibit a den- 
sity jump of about ~ 2, whereas the analytical solution gives 
a compression factor of y + 1/y - 1 -4. These results are 
in accordance w i th previous findings ( Rosswog & Price 1120071 



in accordance w r tn previous findings (Kosswog & FnceluUU/t 
ISpringelll20T0bT: iHefi & SDringell(2010h~ and indicate that for 
3D SPH simulations of the Sedov-Taylor point explosion prob- 
lem the simulation profiles converge to the analytical solution 
as th e resolution is increased , with approximately ~ 345 3 parti- 
cles (iRosswog & Pricel l2007) being required to fully resolve the 
shock front. 

Fig. [9] can also be used to verify the behavior of the indi- 
vidual time-step algorithm. For problems involving ve ry strong 
shocks, as demonstrated by ISaitoh & Makino I d2009l) . individ- 
ual particle time-steps must be properly restricted so to avoid 
that in the proximity of the shock front they do not satisfy the 
local Courant condition, thus leading to inaccuracies in the in- 
tegration. In fact, SPH simulations using individual time-steps, 
but without an appropriate limiter, will fail to predict the e x- 
pected solution profile (see Fig. 3 of lSaitoh & Makiriol d2009l) ). 
The profiles of Fig. [9] demonstrate the correctness of the algo- 
rithm used to update the time-steps, although different from that 
devised bv lSaitoh & Makino I d2009l) . 

The present scheme adopts individual particle time-steps Af,- 
that can vary in pow er-of-two subdivisions of the largest allowed 
time-step A?o > Af, dHernquist &Katz 111989b . At each step, par- 
ticles whose time bin is synchronized with the current time are 
defined as 'active' and their hydrodynamic quantities, as well as 
their smoothing lengths and densities are consistently updated. 
However, non-active particles j, that are neighbors of an active 
particle i are defined as 'low-order active' particles and their hy- 
drodynamic variables, as well as their time-step constraints but 
not their forces, are recalculated. This criterion is applied regard- 
less of the local shock conditions and no particular conditions are 
imposed on Af, which depend on Af,. 

3.4. The blob test 

The blob test is another hydrodynamic test where results of stan- 
dard SPH diff er significantly from th o se produced by grid based 
simulations dAgertz et al. I 120071: iRead. Havfield & Agertz I 



201C : ICha. Inutsuka & Navakshinl 1201 Ot IHefi & Spring 



201C ; iMurante et al. 112011 ). The test consists of a gas cloud of 
radius R c / placed in an external medium ten times hotter and less 
dense than the cloud, so as to satisfy the pressure equilibrium. 
A large enough wind velocity is given to the hot low-density 
medium so that a strong shock wave strikes the cloud. The 
interaction of the cloud with the supersonic medium produces a 
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Fig. 9. Radial density profiles of the 3D Sedov blast wave test 
at f = 0.06. The solid black line indicates the analytic solution, 
while the simulation profiles are obtained by averaging for each 
radial bin SPH densities calculated from the particle distribu- 
tions over a set of grid points located at the surface of a spheri- 
cal shell and uniformly spaced in angular coordinates. All of the 
runs were performed using = 524, 288 particles. 



number of effects that are of interest in an astrophysical context, 
such as gas stripping and fragmentation. 

Initially, the blob is perturbed by the development 
of Richtmyer-Mesh kov and Rayleigh-Taylor instabilities 
dAgertz et al. I l2007h . After wards, large-scale (~ R d ) KH 
instabilities are created at the cloud surface because of the 
shear flows due to the supersonic wind. This non-linear phase 
is supposed to dev elop over a time-scale r er ~ 2R c i -Jxl vw 
dAgertz et al. II2007I) , where x is the density contrast, after which 
cloud disruption will take place. 

In order to investigate the capability of the AC-SPH code 
to properly follow the hydrodynamic of the blob test, we com- 
pare results extracted from a set of SPH simulations realized 
with the same initial conditions but with different numerical 
parameters. The numerical setu p of the test is the same as in 
IRead. Havfield & Agertzld2010l) . to which we refer for more de- 
tails. A spherical cloud of radius R c i = \91kpc is placed in a 
periodic, rectangular box of size \L X , L y , L z ) = {2, 2, %}Mpc. The 
cloud has density p c \ = 4.74 x 10~ 33 grcm~ 3 = xPexi and tem- 
perature T c i = \0 b K = T ex ,lx- The ambient medium has density 
and temperature so that^- =10. The cloud is initially located at 
{1,1,1 }Mpc and the ambient medium is given a wind velocity 
vw - \,000kms~ l , so that for an adiabatic index y — 5/3 its 
Mach number is M — 2.7. An HCP lattice of equal mass parti- 
cle is constructed in order to satisfy the above density require- 
ments so as to use for this ve rsion of the test, a total number of 
N ~ 1.1 X 10 6 particles, as in lHefi & Springel ld2010l) . Finally, a 
velocity perturbation is imposed at the cloud surface in order to 
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trigger the development o f an instability; amplitude and m odes 
are given in appendix B of lRead. Havfield & Agertzl (120101) . 

We compare results from SPH simulations where three dif- 
ferent spline kernels have been used : CS, CRT and M5. We use 
128 neighbors for the simulations with the CS and CRT ker- 
nels and 220 for the run with the M5 kernel, so that the ratio 
77 is the same for all the runs (77 ~ 1.5). For the CS kernel, 
we run a simulation where the AC term (IT~8b is absent in the 
equation of thermal energy evolution (CS-128 NOAC), and three 
other simulations (CS-128, CRT-128, M5-220) in which the AC 
term dT8l is incorporated in the energy evolution equation. In 
the following, the simulation AV parameters are set as follows: 
\a min ,a max ,l d ) = {0.1, 1.5,0.2} 

For several runs Fig.[10]shows the gas density maps at three 
different times: t — 1 , 2, 3; the time being in units of tkh ~ 1 -6T cr 
(lAgertz et al. 1120071) . The maps have been evaluated according 
to the SPH prescription on a 2D grid of (200x800) points in the 
central yz plane located at x — L x /2. The top panels of Fig. \W\ 
are for the standard NOAC run CS-128. It can be seen that 
in this case, unlike what was found in mesh-based simulations 
(lAgertz et al. 112007b . the cloud survives the impact of the super- 
sonic wind striking its surface and there is no disruption. This is 
due to the absence of fluid instabilities developing at the cloud 
surface. If they were present they would in turn lead to the strip- 
ping of material and to the cloud break-up. This code behavior 
is similar to what was found in previous findings (|Agertz et al.l 
120071: iRead. Havfield & Agertz iboiOHHefi & St)ringelll2010h . 

On the contrary, incorporating the AC diffusion term in the 
SPH equations leads to a significant improvement in the code 
capability to properly model the cloud evolution. This can be 
seen from the middle and bottom row panels of Fig. [10] in which 
the density maps are displayed for the CS-128 and M5-220 runs. 
The simulation performed using the CRT kernel is not shown 
since its performances are quite similar to those obtained using 
the CS kernel. 

This behavior is consistent with the expectation that intro- 
ducing AC removes the numerical effects which in SPH pre- 
vent the treatment of contact discontinuities when large density 
jumps are present, and thus the inconsistencies that suppress the 
growth of the instabilities. The AC-SPH formulation presented 
here can therefore, at least qualitatively, correctly follow the 
time evolution of the cloud as in the other SPH schemes which 
have been recently proposed (IRead. Havfield & Agertzl 1201 ot 
Hefi & Springelll2010t fMurante et al. 11201 it ISaitoh & Makino l 
20121) . 

To quantify the code performances in a more quantitative 
way , Fig. QT|(left panel) shows the mass loss of the cloud as 
a functio n of time for the diff erent runs. We follow previous def- 
initions dAgertz et al. 1120071) and a gas particle is defined to be 
still member of the cloud if at the considered epoch its gas and 
temperature fulfill the conditions p > 0.64p c ; and T < 0.9T ext . 
The standard SPH results show a cloud that is not disrupted and 
its mass, after an initial transient, stays constant at about half of 
the original value. A striking result is instead given by the runs 
employing the AC term. For these simulations there is now a 
large mass loss rate occurring at early times, followed by com- 
plete cloud disruption at t Z 3tkh- 

Note that the mass loss rate of the the AC-SPH runs does 
not depend in a significant way on the choice of the kernel. This 
suggests that the cloud disruption is driven by large-scale in- 
stabilities and is relatively insensitive to small-scale perturba- 
tions. Given the similarities displayed in Fig. Q~T] left panel, by 
the mass loss rates of the AC-SPH simulations employing dif- 
ferent kernels, when referring to the left panel we will adopt the 



term mass loss rate to generically indicate the behavior of these 
curves. 

Clearly, in the test considered here the new AC scheme is 
now capable of properly removing the surface effects, which are 
present across the contact discontinuity in the standard SPH ver- 
sion, that artificially suppresses the growth of hydrodynamic in- 
stabilities. However, a comparison of the mass loss rate with the 
corresponding one produced using the new mo ving-mesh code 
AREPO in simulations of equivalent resolution (ISpringel 11201 lL 
Fig. 10), shows that for the AC runs presented here the mass de- 
pletion of the cloud occurs much faster. This discrepancy sug- 
gests that the processes of heat diffusion which in the adopted 
numerical scheme are mediated via the AC parameter a c , should 
be somehow constrained by a physically motivated mechanism 
which has not been previously considered in the discussion of 
Sect. 12.31 This mechanism should be introduced with the pur- 
pose of avoiding a heat transfer mechanism, as governed in the 
code by Eq. ( fT8l , that is overly diffusive. 

In order to evaluate the relative effectiveness of heat and 
momentum transport, in the theory of heat transfer the Prandtl 
number Pr is defined as the ratio between the kinematic vis- 
cosity v and the thermal d iffusion coefficient D: Pr = v/D 
dBlundell & Blundell 1 12006b . For gases, the transport coeffi- 
cients for the transport of heat and momentum are nearlv equal 



and the Prandt l number is of order unity, with Pr 



when 



y = 5/3 dBlundell & Blundell II2006I) . This suggests that a con- 
straint on the particle AC parameter af can be obtained by set- 
ting 



Pr 



vav 
Dac 



cx av v av 
I IJ 



(43) 



where the equivalent kinematic viscosity coefficient vav 
due to AV is given approximately by vf v ~ jQaf v vf^rij 
dLodato &Pricell20lol) . and the numerical heat diffusion coef- 
ficient Df c has been estimated according to Eq. (f2~Tb . Note that 
the reverse of Eq. (l43l is relatively unimportant, since the source 
term S c , is driven by a second order derivative, whilst for AV is 
a first order derivative which determines Si. 

It can be seen that from Eq. (1431 . for the AC signal veloc- 
ity adopted here, the constraint on the numerical heat diffusion 
becomes important in presence of supersonic flows. Moreover, 
the condition Pr ~ (9(1) is valid only for gases, where the trans- 
port of momentum and energy occurs simultaneously. For liq- 
uids, according to the dominant mechanism of heat conduction, 
the value of the Prandtl number can va ry by several ord ers of 
magnitude among different substances dPimotakis 1120051) . The 
condition d43l should, therefore, be considered problem depen- 
dent. 

A rigorous procedure to derive the upper limits on the set 
of parameters Ja?J at any given timestep would require to first 
define for the particle i the numerical kinematic viscosity and 
heat diffusion coefficients 



<vf v > 
<D AC > 



i(£Lj/i v2v < +2V < v - v <) 



(44) 



where the bulk visc osity 4" is 5 /3 of the shear viscosity rj = vp 
dLodato & Price 11201 Oh and it is understood that the SPH expres- 
sions should be used for the operators at the denominators. 

The upper limits on the parameters JaJ'} are then given by 
the conditions 



< Df > < - < v? v 



(45) 
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Fig. 10. Density maps of the blob test in the central plane x - L x /2 at t = 1, 2, 3 for SPH runs with (CS, M5) and without (CS 
NOAC) the AC term. Time is in units of tkh ~ 1 -6T cr and the number accompanying the kernel label indicates the number of 
neighbors of the run. Axis units are in Mpc. 
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Fig. 11. Mass loss of the cloud versus time for the blob test simulations. Particles are defined as cloud member if their temperatures 
and densities fulfill the conditions < Q.9T ext and p,- > 0.64p e /. Left panel: simulations have been performed incorporating the 
AC term into the SPH equations and using the kernels CS, CRT and Ms- For the CS kernel a simulation is also shown as reference 
without AC (NOAC, solid black line). Right panel: for the CS kernel simulations with additional constraints on the particle AC 
parameters a c have been perf ormed usin g N = 10 6 (N6) and N = 10 5 (N5) particles, see text for more details. Black squares have 
been extracted from Fig. 10 of lSpringell (1201 ll) and indicate the behavior of the cloud mass versus time in a similar test performed 
using the moving-mesh code AREPO and with a number of (64 x 64 x 128) resolution elements. 
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which must be simultaneously satisfied by all the parameters. 
In place of the procedure just described, we adopt here a simpler 
approach and using Eq. d43l we estimate the upper limits on jorpj 
as 



of < ^of MAX^vf |/|vf |} , 



(46) 



where the notation j < means that the maximum is taken 
form all the neighbors j that satisfy the condition Vy ■ r, y < 0. 
This is to be consistent with the definition of AV, for which n,y is 
non-zero only for approaching particles. Clearly, this definition 
does not guarantee an accurate estimate of the constraint on the 
af parameter, nonetheless the use of the maximum among all the 
approaching pairs should provide a floor value for the effective 
constraint. 

For the CS kernel, we then performed a simulation identical 
to the one previously discussed, but now imposing the additional 
constraint d46i l on the a c parameters. This simulation is labeled 
as MAX N6 and the corresponding mass loss curve is displayed 
in the right-hand panel of Fig. QT| In order to assess resolution 
effects, we ran a mirror simulation now using N - \0 5 particles 
(MIN N5). Moreover, the uncertainties associated with the use 
of Eq. d46b can be estimated by looking at the mass loss rate 
of a simulation (MIN N6) where instead of Eq. (l46l i the same 
equation was used but the constraint is derived using the MIN 
operator instead of the MAX one. 

For comparative purposes, in the right-hand panel of Fig.fTTI 
the mass loss of the cl oud (black square s), as found in a similar 
blob test performed bv lSpringeTl (1201 ll) using the new moving- 
mesh code AREPO, was also inserted. The number of resolution 
elements is 64 x 64 x 128, so that the resolution is approximately 
equivalent to that of the tests displayed here. 

A striking result seen from the behavior of the mass loss 
rates is that, introducing the constraint (|46| | on the numerical 
heat transfer in the simulations, the mass depletion of the cloud 
is now in better agreement with previous results and in particular 
with the cloud mass evolution produced in the blob test by the 
Voronoi-based code AREPO. 

Uncertainties associated with the procedure adopted to esti- 
mate the constraint on the a c parameters are of limited impact, 
as can be assessed by the differences between the mass loss rates 
exhibited by the run MAX N6 and MIN N6. Moreover, for the 
run MAX N5 the mass loss rate is much stronger than in the 
case MAX N6. This is indicative of numerical diffusion effects 
which dominate the blob evolution. Clearly, a numerical simu- 
lation with N - 10 1 particles should be required to clarify this 
issue; however the agreement with previous results of the mass 
loss for the run MAX N6 suggests that convergence is already 
being achieved using = 10 6 particles. 

Finally, both the runs MAX N6 and MIN N6 exhibit at 
t i 3tkh a remnant cloud mass which is of the order of ~ 
20 - 30% of the initial mass. Although such a result can still 
be due to the use of a constraint on the numerical heat diffu- 
sion which still needs to be refined, we note that such remain- 
ing masses are equally present in SPH blob tests which have 
adopted, with the purpose of avoiding t he problems of standard 
SPH, completely dif ferent approaches dHefi & SpringellfeOlOt 
iMurante et aL~ll201ll) 

We, therefore, suggest that this underestimate in the stripping 
rate of the blob mass is not due to the AC scheme, but depends 
rather on the SPH formulation and is associated with the intrinsic 
errors in gradient estimates. These errors, in turn, lead to the 
suppression of the small-scale instabilities. Such a issue will be 



discussed further in detail in Sect. l4.2l dedicated to the Rayleigh- 
Taylor instability. 

To summarize, SPH simulations which incorporate the AC 
scheme can successfully be used to accurately model the blob 
evolution. However, the results of the tests indicate that it is nec- 
essary to properly constrain the AC parameter a c in order to 
avoid unphysical heat diffusion when strong supersonic flows 
are present. 



4. Gravity tests 

To validate the performances of the AC signal velocity d25l >. sev- 
eral hydrodynamical test problems were investigated in which 
the self-gravity must be necessarily taken into account to prop- 
erly model the system evolution. We first consider the 3D col- 
lapse of a cold gas sphere initially at rest, this test being cus- 
tomarily used in SPH to judge the code capability when strong 
shocks are present, such as those occurring during the formation 
of self-gravitating structures. We then examine a 2D version of 
the classic Rayleigh-Taylor instability test and finally the new 
code is used to assess the behavior of cluster entropies in cos- 
mological structure formation. 



4. 1 . Cold gas sphere 

A standard hydrodynamical test for SPH codes in which gas- 
dynamics is modeled including self-gravity is the 3D col- 
lapse of a cold gas sphe re, also comm only recognized as 
the 'Evrard' collapse test dEvrardlll988l) . The test follows in 
time the adiabatic collapse of a initially c old gas sphere and 



has been widely used by many a uthors (jHer nquist & Katz 
I 1989 I : ISteinmetz & Mueiksrl 119931: IWadsley Stadel & Ouinn 



2004tlSDringell 



iHefi & Springe! 



20051: IWetzstein et al. O2009t ispringel 11201 Obt 
120101 VI 1) as a standard test for SPH codes. 



The gas cloud is spherically symmetric and initially at rest, 
with mass M, radius R, and density profile 



p(r) = 



M 1 



2kR 2 r 



(47) 



The gas obeys the ideal gas equation of state with y — 5/3 
and the thermal energy per unit mass is initially set to u — 
Q.05GM/R. The SPH simulations are performed using units for 
which G = M — R — 1 and th e chosen time unit is the cloud 
free-fall timescale t /f = (tt 2 /8) ^R 3 /GM = tt 2 /8. 

With these initial conditions, the pressure support of the gas 
sphere is negligible and the cloud begins to collapse until a 
bounce occurs in the core with a subsequent shock wave propa- 
gating outward. Most of the kinetic energy is converted into heat 
at the epoch of maximum compression of the gas, which occurs 
at t ~ 1.1. The initial conditions setup is realized by stretch- 
ing the radial coordinates of a glass-like uniform distribution of 
= 88, 000 equal mass particles located within a sphere of unit 
radius, so as to generate the density profile of Eq. (|47] >. 

To construct the set of SPH simulations, four tests have been 
run using the kernels CS and M5. The number of neighbors is 
set to 100 for the M5 runs and 50 for the CS ones. To assess 
convergence properties in the radial profiles of the considered 
hydrodynamic variables, the M5 run has been replicated using 
200, 000 particles (M5 - 200k). All of these runs were performed 
incorporating in the SPH equations the AC term (fT8l , for the CS 
kernel a reference run with the AC term disabled has been con- 
sidered (CS-NOAC). The gravitational softening length is taken 
as Eg = 0.02. 
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Fig. 12. Radially averaged profiles at t — 0.8 of the Evrard collapse test. Clockwise from the top left panel: profiles of density, 
entropy, artificial conductivity parameter a c and radial velocity. Curves with different line styles and colors refer to SPH runs 
performed using different kernels and with (AC) or without (NOAC) the AC term in the SPH energy equation. The black solid lines 
indicate the ID PPM reference solution of lSteinmetz & Mueller] (Il993l) . 



For each test case, we perform simulations with different 
settings of the AV parameter lj. However, introducing a time- 
dependent AV scheme in SPH reduces numerical viscosity ef- 
fects and in particular, for the test investigated here, produces a 
radial entropy profile at the shock front in better agreement with 
the ID PPM reference solution. For this reason, the results for 
different test cases are presented for those runs where a standard 
AV formulation has been used ( runs AVq of VI 1). This is done in 
order to highlight differences in the simulation results due to the 
use of different kernels and of incorporating the AC term in the 
SPH thermal energy equation. This choice for the simulation pa- 
rameters allows one to compare the test calculations performed 
here with previous results presented in Sect. 5.2 of VI 1 . 

The average radial profiles at t — 0.8 of density, entropy, 
time-dependent AC parameter and radial velocity are displayed 
in Fig. [12] The black solid lines indicate th e profi les of the 
ID PPM calculation of lSteinmetz & Mueller I d!993l) . Broadly 
speaking, one expects to see some differences between the pro- 
files of the NOAC simulation and those of the other runs. This 
should be valid in particular for the entropy profiles. However, 
at the considered epoch, Fig. [12] shows that the entropy profile 



of the NOAC run is still very similar to the others. This occurs 
because at t — 0.8 the shock front propagation is still in the early 
phase and the smoothing in the thermal energy due the AC term 
is not yet significant. 

The behavior of the a c profiles is in accordance with what 
expected. In particular, the profiles exhibit a peak in correspon- 
dence of the shock front location, which occurs approximately 
at r ~ 0.18, and a decay as one moves away from it. There is 
a dependence on the kernel shape and a weak dependence on 
resolution. The former is interpreted as a consequence of an im- 
proved mixing capability due to the increased order of the ker- 
nel. Note that there is a rising in the profiles as one moves to- 
ward the sphere center, this occurs because of the presence of 
a temperature gradient. However, in Eq. ( [T8l this radial depen- 
dence is of no effect since behind the shock front the sphere 
quickly achieves a hydrostatic equilibrium and vrf{grav) 
forr < 0.1. 

The epoch examined in Fig.[T2lhas been chosen so as to com- 
pare the simulation profiles w i th the c orresponding ones pre- 
sented in Fig. 40 of ISpringell (1201 Obi) and realized using the 
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Fig. 13. As in Fig.Q~2l a closer view at t — 0.9 of the radial profiles of density and entropy in the proximity of the shock front. 



new moving-mesh code AREPO. In particular, a visual inspec- 
tion shows that ahead of the shock front the accuracy of the en- 
tropy profiles shown in Fig. [12] can be considered intermediate 
between the profile produced using the AREPO code and the one 
re alized with the SP H code Gadget-2 ( right panels in Fig. 40 
of lSpringell(l2010bl) . test cases 'moving-mesh' and 'SPH', re- 
spectively). Note, however, that a strict comparison between the 
different simulation profiles is di fficult because the number of 
resolution elements employed by ISpringell d2010bl) in the cold 
gas sphere tests are lower by about ~ 1/3 than the one used here. 
In particular, the smaller shock broadening found here can be as- 
cribed to a resolution effect. 

Finally, Fig. [13] shows a closer view of the radial density 
and entropy profiles in the proximity of the shock front at t — 
0.9. The NOAC simulation exhibits now the worst accord with 
the PPM reference solution, while for the high-resolution run 
M5 - 200fc the improvement in accuracy is minimal with respect 
to the other AC runs. This suggests that to obtain a significant 
improvement in the adherence of the simulation profiles to the 
PPM solution one mu st use a much high er number of simulation 
elements ( say > 10 6 , ISpringeil d2010bl) ). 

A noteworthy feature of the profiles of Fig. [13] is that now 
the AC simulation profiles are in much better agreement with the 
PPM profiles. This is interpreted as a numerical effect in which 
the pre-shock entropy, which is generated by the AV implemen- 
tation ahead of the shock front owing to the converging flow, is 
now strongly reduced because the presence of the AC term in 
the SPH energy equation now removes this excess of internal 
energy. This in turn implies, behind the shock front, a thermal 
energy profile closer to the PPM solution profile. 

Comparisons of the profiles of Fig. [13] with the correspond- 
ing ones displayed in Fig. 4 of VI 1 show that the dependency of 
the simulation profiles on the adopted AC scheme are stronger 
than the ones due to the time-dependent AV formulation. In fact, 
these results are almost unaffected if one performs SPH simu- 
lations using now a time-dependent AV scheme in place of the 
standard one. 

Finally, two runs have been performed using the CRT and 
LIQ kernels and a number of neighbors set to N sp i, = 50. The 
profiles of the two simulations have not been shown here to avoid 



overcrowding in the panels, however the performances of the 
CRT run are quite similar to those of the CS one, whilst the pro- 
files of the LIQ simulation exhibits very poor consistency prop- 
erties with respect to the reference PPM solution profiles. 

4.2. 2D Rayleigh-Taylor 

The Rayleigh-Taylor (RT) instability a rises when a heavier fl uid 
is placed on the top of a lighter fluid (IChandrasekhar|[l961l) in 
presence of an external gravitational field. The fluids are in pres- 
sure equilibrium against the external field and in this configu- 
ration the system is unstable in the presence of small perturba- 
tions at the interface, the lighter fluid will then begin to rise and 
the denser one to fall. This process leads to the development 
of characteristic finger-like structures before the fluids enter the 
non-linear phase where they mix together completely. In order 
to evaluate the ability of the AC-SPH code to properly describe 
the evolution of RT instabilities, we consider a 2D version of the 
test. The initial conditio ns are chosen to be similar to the numer- 
ical test implemented bv lAbell I d201 ll) to validate its new rpSPH 
formulation of SPH equations, so as to consistently compare the 
results. The computational domain consists in a 2D box with co- 
ordinates x e {0, 1 /2), y e {0, 1). The boundaries are periodic in 
x and reflecting in y. The density is p\ -2 at the top and pi-\ 
at the bottom, with a density profile 



p(y) =P2 + 



(Pi - Pi) 



[l. + exp{-2(j-0.5)/A v }] 



(48) 



where A v = 0.05. This ensures a smoothing in density at 
the interface y - 0.5 that a llows a consistent numerical behav- 
ior dRobertson et al. ||2010|) . This density profile is realized by 
constructing a HCP lattice of N = 620 2 equal mass particles 
in which the spacing is varied according to the procedures de- 
scribed in Sect. 13. l.TI until the profile d48l i is satisfied. 

The pressure at the interface is set to Pq = pi/y = 10/7, 
where y = 1.4 and it varies with y according to P(y) = Pq - 
gp(y)(y- l/2),g = 1/2 being the external acceleration, so that the 
system is initially in hydrostatic equilibrium. For the particles 



21 



R. Valdarnini: Hydrodynamic of an SPH code incorporating an artificial conductivity term 



that satisfy the condition \y - 0.5 1 < 0.2 a velocity perturbation 
is applied in the y direction given by 



v y (x,y) = 6v y {l +cos[8tt(x+ 1/4)]} 

{l+cos[57r(j-l/2)]}/4, 



(49) 



where 6v y =0.1. 

We present results for the CS, CRT, Ms and additionally, 
for reasons which will be discussed later, we consider also the 
Me kernel dPrice II20T2I sect. 2.3). The number of neighbors for 
these runs is set to N sp h = 32,32,50, 110, respectively. As in 
the other test cases considered here, for comparative purposes, 
we ran a simulation (CS-NOAC) in which the AC term in the 
thermal energy equation has been disabled. The following set- 
tings have been used for the AV parameters: {a m i n , <Xmax> ld\ = 
{0.1,1.5,0.2}. 

In Fig. [14j at the time t = 0.3, we present the 2D density 
maps of the different RT tests. As expected, standard SPH is un- 
able to correctly capture the development of the RT instability, 
as indicated by the leftmost panel of Fig.[14](CS kernel, NOAC 
run). However, from this viewpoint the improvement is mini- 
mal even for the corresponding AC version of the considered run 
(CS-32). On the contrary, a significant improvement is obtained 
if one uses the M5 kernel or th e next in order M^, which is a 
quintic polynomial dPrice II20T2I) . In fact, in the rightmost panel 
of Fig.[14](M6 run), the center of mass of the RT spikes appears 
located at a lower position than in the M5 run and therefore the 
convergence might still not be achieved even for the run. 

This strong dependence of the code performances on the ker- 
nel order demonstrates that, for the RT test, accuracy in pres- 
sure forces is a fundamental issue. These results are consistent 
with those of sect. l3.1l and indicate that the poor performances of 
SPH when handling KH or RT instabilities, in particular for very 
subsonic flows, are mainly due to the leading errors in the mo- 
mentum equation. These errors are reduced as the order of the 
kernel is increased, therefore implying an improved accuracy in 
pressure force estimates and thus a lower velocity noise. This, 
in turn, implies a better capability to capture the growth of the 
instability. 

Similar results h a ve r ecently been obtained by 
iMcNallv.Lvra&Passvl (120111) . who ran a suite of 2D KH 
simulations using carefully crafted initial conditions with the 
goal of assessing different hydrodynamic code capabilities. 
The authors conclude that for SPH the code performances are 
strongly related to the order of the kernel employed in the 
simulation and thus to the accuracy of the gradient estimates. 

Moreover, the si mulation behavior of Fig. [T4l is consisten t 
with the findings of iGarcfa-Senz. Cabezon. & Escartfn I (120121) . 
The authors present a new formulation of SPH in which gradi- 
ents are estimated using a tensorial approach which conserves 
both linear and angular momentum. The authors demonstrate 
that using this higher order gradient estimator the interpolation 
of physical quantities is significantly improved with respect to 
the standard method. In particular, several tests showed that the 
code is now able to successfully capture the development of KH 
and RT instabilities. 

To summarize, the results of this section demon- 
strate that SPH performances to model hydrodynamic 
subsonic instabilities depend critically on the accuracy 
of the gradient estimator. The f ormula tion proposed by 
IGarcfa-Senz. Cabezon. & Escartfn I (1201 2l) looks particu- 
larly promising in this aspect to overcome the present SPH 
difficulties. 



Finally, it must be stressed that the simulation results of 
Fig. [14] have been obtained setting 6v y = 0.1. If the same 
test had been performed using 6v y = 6.01, as in the R T test 
of the new moving-mesh code AREPO dSpringel 11201 Obi Sect. 
8.8), the velocity noise would have suppressed the growth of 
the instability even for the M(, run . Note, that in their RT test, 
IGarcfa-Senz. Cabezon. & Escartfn I (12012b were able to recover 
the growth of the instability using the same amplitude of the ve- 
locity perturbation, 8v y = 0.01. 

4.3. Cluster comparison 

In this section, a set of cosmological cluster simulations is con- 
structed using the standard SPH formulation as well as the 
new AC-SPH version. We only consider non-radiative, or 'adia- 
batic' simulations, in which the hydrodynamic is modeled ac- 
cording to t he formulation presented in Sect. El Previous in- 
vestigations dFrenk et al. ll999HWadslev Stadel & Ouinn 12004 



Wadslev. Veeravalh & Couchman I I2008T iMitchell et al 



2009; 



Springell l2010bir showed that for the same initial conditions 



there are systematic differences between results extracted from 
cluster simulations produced using standard SPH and Eulerian 
AMR codes. Specifically, the level of central entropy is found 
to be lower in SPH simulations than in the corresponding AMR 
runs. The latter simulations are characterized by the presence of 
a flat entropy core, whereas the radial entropy profile produced 
by SPH runs increases steadily with radius. As already outlined 
in the Introduction, the origi n of this discrepancy has been the 
subject of an intense debate dAgertz et al. 1120071: IMitchell et al."1 
120091: iRead. Havfield & Agertz 112010b which has given as main 
explanation the different degrees of numerical mixing present in 
the two codes. 

To assess the capability of the AC-SPH code in solving this 
issue, we then run several cluster simulations. A comprehensive 
study of the differences between the hydrodynamic variables of 
cluster simulated using the standard and the AC version of the 
SPH code will be presented in a forthcoming paper. Here, we just 
show the final radial entropy profile of the simulated clusters as 
the main variable that can be used to test the effectiveness of the 
AC-SPH formulation for solving the entropy problem in cluster 
cores. 

A detailed description of the procedures used to construct the 
simulated cluster sample is given in sect. 2 of VI 1, we provide 
a brief summary here. The simulations were carried out assum- 
ing a spatially flat ACDM model, with matter density parameter 
£l m = 0.3, vacuum energy density Qa = 0.7, baryonic density 
Q b = 0.0486 and h = H/lOOkmsec^Mpc- 1 = 0.7. The scale- 
invariant power spectrum is normalized to erg = 0.9 on an 8 h~ l 
Mpc scale at the present age to. An N-body cosmological simu- 
lation involving only gravity was first run with a comoving box 
of size L2 = 200/i~'Mpc starting from an initial time f,„ down to 
the final epoch to- At this epoch clusters of galaxies are identi- 
fied in the simulation box as groups of particles that are associ- 
ated with overdensities approximately in excess of ~ 200O m 6 . 
Several of these clusters are then resimulated individually using 
the AC-SPH code here described, coupled with a treecode grav- 
ity solver. 

Initial conditions for a specified cluster are generated as fol- 
lows: a spherical region with origin located at the cluster cen- 
ter is populated with a high-resolution (HR) grid, the radius of 
the HR sphere is such that it contains all of the original parti- 
cles identified at / = 1 as cluster members. A gas and a dark 
matter particle is associated to each grid node, whose positions 
are perturbed according to the random realization of the original 
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h CS-32-N0AC 




Fig. 14. Density maps of the 2D Rayleigh-Taylor tests are shown at time t — 3 for standard SPH (NOAC) and AC-SPH simulations 
performed using N = 620 2 parti cles with diffe rent kernels and neighbor numbers. The density distributions can be compared directly 
with the map shown in Fig. 4 of lAbell I (1201 ll) . top right panel, as the adopted initial conditions for the tests are the same. 



cosmological simulations and only those particles which are lo- 
cated within the HR sphere are kept for the hydro run. The HR 
particles are surrounded by a low-resolution shell of dark matter 
particles, extracted from a grid with spacing twice that of the HR 
grid, the shell is introduced with the purpose of mimicking the 
effects of tidal forces. The grid spacing is chosen such that at the 
end of the procedure a cluster is simulated with N gas ~ 220, 000 
gas particles and Ndm ~ N gas dark matter particles in the inner 
HR sphere, whereas N%% ~ Ndm are use d m the low-resolution 
shell. Particle masses are assigned according to the values of O m 
and fib- 

The clusters selected to be re-simulated hydrodynamically 
were chosen with the following criterion. Originally (VI 1), the 
procedure previously described was repeated two more times in 
order to construct two new cluster samples extracted from cos- 
mological simulations with box sizes L4 = 2L 2 and L% = 2L4. 
The three cluster samples were then combined to construct a 
final sample S a u of ~ 160 clusters covering nearly a decade 
in cluster masses. Cluster members of sample Sail were then 
ranked according to their dynamical state, as measured by an 
appropriate statistical indicator. Two sub-samples of four test 
clusters each, denoted respectively by Q and P, were then con- 
structed by extracting from sample S a ii those clusters with mem- 
bership criterion for sample Q (P) of being in a fully relaxed 
(perturbed) state. These two sub-samples were then used to in- 
vestigate the effects of AV in hydrodynamic SPH simulations of 
galaxy clusters (VI 1). Here the four test clusters of sub-sample 
Q were chosen to investigate the performances of the AC-SPH 
code. This choice being motivated by the need to proper disen- 
tangle the amount of entropy mixing produced during merger 
processes from that specific to the numerical method in the fi- 
nal cluster entropy. The clusters of sub-sample Q were carefully 
selected on the basis of their dynamical state and are among the 
most relaxed of sample S a ii- Therefore, for these clusters, differ- 
ences in the final radial entropy profile between standard and AC 
runs can be safely ascribed to the new AC term implemented in 



the code. The AV simulation parameters of the simulations are 
{a min ,a max ,l d } = {0.1, 1.5,0.2}. 

The final entropy profiles for the four test clusters are dis- 
played in Fig. [15] where the cluster entropies have been normal- 
ized to 



5 2 oo - \ 
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15 f b H 



2/3 



(50) 



where ft, = Qb/Q m is the global baryon fraction and Ma = 
(4n/3) A p c r\ denotes the mass contained in a sphere of radius 
r& with mean density A times the critical density p c = 3H 2 /8jtG. 
As is customary, the cluster mass is defined setting A = 200 and 
r2oo denotes the cluster radius. 

The results of Fig. [TBI indicate that, although there is a cer- 
tain degree of scatter between the entropy profiles of individual 
clusters, nonetheless all of the AC-SPH runs consistently exhibit 
much shallower entropy profiles and higher core entropies than 
in the standard SPH runs at r/r>oo ~ 0.1. For a given cluster, the 
difference in the levels of central entropies produced by the two 
codes is about of a factor ~ 4, and now the values of central en- 
tropies are comparable with those pr oduced using AMR codes 
in cosmological cluster simulations (IVoit. Kav & Brvanl 20051) 
and in idealized cluster binary mergers dMitchell et al 112009b . It 
should, however, be stressed that in Eulerian hydrodynamics it 
is the numerical scheme that forces the fluid to be mixed below 
the minimum cell size. This suggests that AMR simulations of 
galxy clusters might overesti mate the correct l evel of core en- 
tropy because of fluid mixing (ISDringelll20i0bh . To summarize, 
the development of entropy cores in the AC-SPH runs is clearly 
a consequence of the heat diffusion term (fT8T l. now present in the 
energy equation. This term acts to redistribute the internal en- 
ergy produced in the shocks during dissipative processes, so that 
the subsequent entropy mixing leads, in the end, to higher levels 
of core entropies. 

Similar results were obtained by 

IWadslev. Veeravalli & Couchman I (120081) . who similarly 
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Fig. 15. Final radial entropy profiles as a function of r/r2oo 
for the four relaxed test clusters. The gas entropy S(r) = 
kBT(r)/fxm p pV 3 is plotted in units of 5200 Different line styles 
refer to different clusters, the numbers indicating the cluster 
membership in the cosmological ensemble from which they have 
been extracted. Thin (black) lines are for the profiles of runs per- 
formed using the standard SPH formulation, while thick (red) 
lines refer to the AC-SPH runs incorporating the AC term. 



added a heat diffusion term to the SPH energy equation, but 
using the prescription (f26b in stead of the Eg. (fl8b a dopted here. 
According to I Wadslev. Veeravalli & Couchman I ([2008), the 
value of the coefficient C tends to be problem dependent and 
for the cluster simulations the best agreement with the entropy 
profiles extracted from the mesh code runs is obtained setting 
C =i 0.1 - 1. Here, the time-dependent formulation (l22l presents 
the advantage of minimizing thermal diffusion away from 
shocks, as demonstrated in the cold gas sphere test by Fig. Q~2l 
where the radial profile of the AC parameter a AC is peaked at 
the shock location. 

The simulations presented in this section have been real- 
ized using adiabatic gas physics. However, a realistic modeling 
of the intra-cluster medium physics must incorporate the possi- 
bility that the gas cools radiatively. Previous simulations (VI 1) 
showed that the presence of cooling leads to the development in 
the inner cluster regions of dense compact cool gas cores and 
subsequently of high levels of turbulence, the latter being pro- 
duced by the hydrodynamical instabilities generated by the in- 
teraction of the compact cool cores with the ambient medium. 

How this scenario is affected by incorporating into the SPH 
equations a numerical heat diffusion term is an issue which can 
be properly addressed only by resorting to numerical simula- 
tions. Here, we outline that two competing effects are expected 
to influence the final level of turbulence in a cluster core: the 
improved capability of the code to describe the development of 
hydrodynamical instabilities, and thus an increase in the amount 
of turbulence, and the reduced availability of cool gas because 
of the presence of material of higher entropy which has been 



brought in the inner core due to the enhanced fluid-mixing prop- 
erties of the code. 

Finally, a resolution study was already performed in VI 1, in- 
dicating for the simulation resolution employed here a substan- 
tial stability in the entropy profiles of the simulations. Although 
the simulations of VI 1 do not incorporate the AC term, we do 
not expect strong variations in the resolution dependence of the 
profiles, given the similarity of the level of core entropies with 
those found using mesh-based codes. 



5. Summary and conclusions 

In this paper, we have presented an SPH numerical scheme 
which incorporates an artificial conductivity term and uses an ap- 
propriate signal velocity for simulation s including gr avity. The 
AC formulation has been introduced by iPrice I ((2008) as a solu- 
tion to the problems encountered by standard SPH to correctly 
follow the development of KH instabilities, due to the incon- 
sistencies of the standard formulation in the description of den- 
sity at contact discontinuity. Here, a suite of hydrodynamic test 
problems is investigated with the purpose of validating the new 
AC-SPH code and to assess its performances when using the 
specifically adopted signal velocity. 

The results of the KH instability test are presented in Sect. 
13.11 in which the code capabilities have been tested by consid- 
ering SPH simulations of the KH test performed using a large 
variety of initial condition set-up and SPH kernels. The se t of 
initial conditions has been chosen as in IValckeetal. I (120101) . so 
as to consistently compare the r esults. These a re in a ccordance 
with the corresponding ones of I Valcke et al. and indi- 

cate that, for the version of the KH test analyzed here, the AC 
implementation is important for the long-term behavior of the 
simulation. 

M oreover, the results of the KH test indicate (fValcke et al 
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2011) that the poor performances of standard SPH to properly 



treat KH instabilities can be explained in terms of two distinct 
effects. The first is a general problem of consistency, which 
for the problem under consideration requires that smooth in- 
terfaces should be present at contact discontinuities, in order to 
obtain numerical convergence. The second effect which in stan- 
dard SPH suppresses the growth of KH instabilities is the lead- 
ing error in the momentum equation, due to incomplete kernel 
sampling dRead. Havfield & Agertz ||201 Oh and quantified by the 
norm E° defined by Eq. d36b . 

To circumvent this problem several proposals have b een 
made dValcke et al. I 1201 Ot IRead. Havfield & Agertz 1 1201 Oh in 
which the standard SPH cubic spline M4 kernel is replaced by 
the new LIQ or CRT kernels with steeper central profiles. These 
kernels present the advantage of being stable against particle 
clumping, so that the number of neighbors can be safely in- 
creased in order to reduce the error E°. 

In this paper, we consider the possibility of reducing sam- 
pling errors by considering higher-order B-spline kernels, 
specifically the M5 or quartic spline kernel. A striking result of 
the KH runs presented in sect. 13.11 is that, for a given value of 
the ratio 77 between the smoothing length and the mean interpar- 
ticle spacing, simulations of the KH test performed using the M5 
kernel have amplitudes of the E° error substantially smaller and 
in line with the behavior of the same quantity for simulations 
employing the LIQ or CRT kernels. A linear stability analysis 
reveals that this result follows owing to the very good stability 
properties of the M5 kernel, in fact the analysis suggests that the 
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clumping instability is absent for values of 77 up to 77 i 2, which 
in 3D corresponds to N sp h ~ 520 neighbors. 

These findi n gs ar e consistent with the recent results of 
iDehnen & Alv I d2012l) who proposed to adopt, wit h the spe- 
cific p urpose of avoiding the clumping instability, the IWedlandl 
111 995ft functions as a new class of kernels. The authors showed 
that in terms of stability and accuracy properties, the quartic 
spline performs extremely well when compared with the pro- 
posed Wedland functions. These results are strictly connected to 
the properties of the kernel Fourier transform, which according 
to the authors must be non-negative to avoid particle clumping, 
and are consistent with the findings of sect. 13.1.31 since increas- 
ing the kernel order both the Z?-spline and the Wedland func- 
tions approach the Gaussian. 

We therefore propose, as a compromise between the need of 
reducing sampling errors while keeping the computational cost 
at a minimum, the use of the M5 kernel with neighbor number 
in the range N sp h ~ 60 - 1 20 as the standard combination which 
guarantees sufficient accuracy in many SPH simulations of as- 
trophysical problems. Moreover, the results of the Sod shock 
tube test of sect. 13.21 demonstrate that in order to obtain simu- 
lation profiles in accordance with the analytic solution, for sim- 
ulations employing kernels with a modified shape the use of a 
much larger number of neighbors than in the case of the M5 runs 
is necessary. 

The results of the gravity tests show that the adoption of the 
AC-SPH scheme significantly reduces, at the level tested in this 
paper, the differences seen in the hydrodynamics between stan- 
dard SPH and grid-based simulations of self-gravitating struc- 
tures. For the cold gas sphere, the entropy profile is in better 
agreement with the PPM reference solution and for the cos- 
mological cluster simulations, a key result is the final level of 
core entropies, which are consistent with those of the central en- 
tropies produced using AMR codes. Thus, it appears that in hy- 
drodynamic simulations where self-gravity is important the AC 
term, accompanied by the proposed signal velocity, plays a key 
role as a mechanism of redistributing thermal energy and hence 
as a source of entropy mixing. 

To summarize, results extracted from simulations of hydro- 
dynamic tests where self-gravity dominates are in much better 
agreement with the corresponding ones obtained using mesh- 
based codes. The results then demonstrate the capability of the 
implemented AC-SPH scheme to properly follow the formation 
of cosmic structures. 

It is worth noting t hat the artifici al heat conduction term was 
originally proposed bv lPrice I (120081) with the purpose of avoid- 
ing the inconsistencies encountered by standard SPH in presence 
of density steps at contact discontinuities. A complementary 
view h as been proposed by IWadslev. Veeravalli & Couchman I 
(12008ft . who introduced the same term, albeit in a different nu- 
merical formulation, with the aim of modeling the level of dif- 
fusion due to turbulence. The two interpretations are not mu- 
tually inconsistent, however the results of the self-gravity tests 
presented in this paper support the view of a heat diffusion term 
which in SPH is capable of mimicking the diffusion due to tur- 
bulence. In a similar fashion, IVioleau & Issa I d2007l) presented 
an SPH scheme to model a free-surface incompressible flow 
which , in ana l ogy w ith 3D Large Eddy Simulations, assumes a 
ISmagorinskv I (1 19631) model for the filtered Navier-Stokes equa- 
tions. 

The results of the blob test simulations demonstrate that the 
instabilities leading to the expected cloud disruption can develop 
only when the SPH energy equation incorporates the AC term. A 
particularly interesting result is that an appropriate limiting con- 



dition must be implemented on the AC coefficients a c in order 
to avoid an unphysical amount of heat diffusion, which in turn 
leads to a cloud disruption which occurs too early. This limiter 
has been identified as given by the Prandtl number and, for the 
AC signal velocity adopted here, it severely limits the amplitude 
of the AC coefficients in the regime of strong supersonic flows. 

AC-SPH simulations of the blob test incorporating now the 
new constraint support this view, since Fig.[TT]shows for the new 
runs a cloud mass-loss rate which is in better agreement with the 
rates obtained from simulatio ns realized usin g a completely in- 
dependent numerical scheme (ISpringel 11201 II) . However, it must 
be stressed that the condition d43l has been calculated for a per- 
fect monoatomic gas with y = 5/3, so that a physically moti- 
vated constraint on the artificial heat diffusion of the simulated 
medium should be considered problem dependent. 

However, the code has still several problems which render 
its use problematic if the development of hydrodynamic insta- 
bilities need to be followed in the regime of subsonic flows. 
The results of the simulations indicate that these shortcomings 
are not due to the AC implementation, but rather are intrin- 
sic in the standard formulation with which gradients are cal- 
culated in SPH and the related errors are subsequently intro- 
duced in the momentum equation. In particular, simulations of 
the RT test show that increasing the kernel order alleviates the 
problems but does not solve t hem. These results are in ac- 
cordance with recent findings dMcNallv. Lyra & Passvl 1201 U 
iGarcfa-Senz. Cabezon. & Escartfn 20121) and clearly demon- 
strate that for very subsonic flows, the poor performances of SPH 
to model hydrodynamic instabilities are strictly connected to the 
code accuracy in gradient es timates. 

Th e formulation of IGarcfa-Senz. Cabezon. & Escartfn I 
d2012l) . which has been proposed with the aim of calculating 
SPH gradients with an as high as possible accuracy while 
keeping the benefits of a Lagrangian formulation, looks in this 
aspect very promising and suggests further investigations along 
the numerical approach proposed in this paper. These would 
be particularly re l evant in the light of the recent results of 
iBauer & Springell d2012l) . who claim that standard SPH fails 
to properly model the regime of subsonic turbulence. They 
reach their conclusions by comparing results extracted from 
simulations of driven subsonic turbulence realized using the 
moving-mesh code AR EPO and GAD GET SPH. Their results 
have been criticized by iPrice I (1201 2l) . for whom the use of a 
time-dependent AV is critical in SPH simulations of subsonic 
turbulence. 

A re-analysis of the IBauer & Springell d2012l) simulations 
using the AC-SPH code presented here, augmented with im- 
proved gradient operators, would then be of fundamental im- 
portance to achieve a deeper understanding of the capability 
of different numerical methods to model subsonic turbulence. 
The latter being expected to have a significant impact in shaping 
the thermodynamic properties of baryons in cosmological haloes 
and, subsequently, the process of galaxy formation. 
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